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1 .  INTRODUCTION 


(1.1)  Background 

The  topics  of  dynamic  loading  of  gear  teeth  and  the 
deflections  of  gear  teeth  due  to  dynamic  loads  have  been 
treated  extensively. 

One  such  work,  presented  by  Cornell  and  Westervelt  [1], 
utilizes  an  improved  version  of  a  model  developed  by  Richardson 
[4].  The  model  generates  the  dynamic  loads  for  a  meshing  gear 
using  a  cantilever  beam  with  a  cam  moving  along  it,  simulating 
the  engagement  and  disengagement  of  the  adjacent  tooth  (see 
Figure  1-1).  These  dynamic  loads  are  then  used  in  a  dynamic 
model  of  meshing  gear  teeth  where  the  two  gear  hubs  act  as 
rigid  inertia  and  the  teeth  as  variable  stiffness  springs  as 
shown  in  Figure  1-2.  Of  significant  importance  in  this 
investigation  is  the  claim  made  by  the  authors  that  the  effect 
of  variable  tooth  stiffness  is  small,  changing  the  dynamic  load 
response  slightly  compared  to  a  system  with  constant  tooth 
stiffness. 

Another  dynamic  load  response  algorithm  was  developed  by 
Wang  and  Cheng  [2-3],  where  they  reported  that  both  the  dynamic 
load  and  the  induced  dynamic  response  are  highly  dependent  on 
the  speed  of  the  moving  load.  In  slow  speed  regions,  the  dyna¬ 
mic  load  response  is  composed  of  a  static  response  which  varies 
with  the  stiffness  of  the  tooth.  Superimposed  on  the  static  is 
an  oscillatory  response  caused  by  the  excitation  of  the  system 
at  the  resonant  frequency.  Wang  states  that  as  the  speed  of 


the  moving  load  increases  and  the  resonant  frequency  of  the 
system  is  approached/  the  dynamic  load  response  becomes  so 
abrupt  that  tooth  separation  occurs.  A  much  smoother  response 
is  generated  when  the  speed  of  the  moving  load  is  increased  and 
becomes  out  of  phase  with  the  system  resonance.  Here  the 
peak  response  is  reduced  significantly  and  actually  becomes 
less  than  that  for  a  static  load.  Examples  of  the  dynamic  load 
variation  obtained  by  Wang  and  Cheng  are  included  in  Figure 
1-3. 

Kasuba  [5]  presents  an  algorithm  which  analyzes  spur 
gearing  under  static  and  dynamic  loading  conditions.  In  his 
analysis,  the  stiffness  of  the  teeth  are  determined  by  solving 
the  statically  indeterminant  problem  of  multi-pair  contacts, 
changes  in  contact  ratio,  and  meshing  gear  deflections.  in 
general,  Kasuba  states  that  to  decrease  the  dynamic  load 
response,  increased  damping  and/or  contact  ratio  can  be  used. 
He  also  noted  that,  in  a  general  sense,  high  contact  ratio 
gears  have  lower  dynamic  loads  than  low  contact  ratio  gears. 

Up  to  now,  the  discussion  of  models  developed  to  determine 
the  dynamic  response  of  gear  teeth  has  been  limited  to  theore¬ 
tical  cases.  Wallace  [9]  in  his  investigation,  uses  finite 
element  analysis  in  conjunction  with  experimental  techniques  to 
study  the  deflection  of  gear  teeth.  He  subjects  a  single 
tooth  to  both  Hertzian  impact  and  general  dynamic  loads,  hoping 
to  define  a  procedure  for  predicting  deformation  distributions 
due  to  dynamic  loads.  The  finite  element  method  shows  good 
correlation  with  experimental  results  obtained  using  a  short 
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(a)  Dynamic  Load  Variation,  r0  «  i,Nt  -  26,  0-  .  e.D  «  .25* 
crn(.1  ln.)(t)n  rn  1S37RPM .(b)n  •  7640  RPM,  and  fc Jo  =  16370 RPM 


DYNAMIC  LOAD  VARIATION 


(b)  Dynamic  Load  Variation,  re  ■  i,  Hr  »  26, 0,  -  6, 0  *  .264 
cm(.1  In.) (a) n  -  1637RPM,(b)n  -  7640  RPM,  and(cjn •  16270 RPM 


DYNAMIC  LOAD  VARIATION 


Figure  1-3 


cantilever  beam  subjected  to  impact  loadings  at  different  posi¬ 
tions. 


Another  important  contribution  to  the  subject  of  gear 
dynamics  was  made  by  Attia  [6).  He  studied  the  effects  of 
including  the  rim  when  performing  a  static  analysis  to  deter¬ 
mine  tooth  spring  constants.  He  concluded  that  the  stiffness 
of  teeth  with  the  rim  included  is  significantly  less  compared 
to  a  variable  cross-section  cantilever  beam  rigidly  fixed  to 
the  gear  body.  With  the  added  flexibility,  the  initial  con¬ 
ditions  of  two  meshing  teeth  are  highly  dependent  on  the 
deflections  of  the  two  previously  engaging  teeth.  This  fact  is 
very  significant,  as  it  will  definitely  affect  the  type  of  load 
experienced  by  the  upcoming  gear  pair. 

Many  of  the  theoretical  models  used  to  predict  the  deflec¬ 
tions  of  gear  teeth,  such  as  those  presented  by  Cornell  [1], 
Wang  [2-3]  and  Kasuba  [5],  make  use  of  tooth  stiffness 
variations  obtained  from  a  static  deflection  analysis.  The 
equations  of  motion  are  expressed  as  functions  of  the  load 
position  only. 

Nagaya  and  Uematsu  [7]  state  that  because  the  contact 
point  moves  along  the  involute,  the  dynamic  load  response 
should  be  considered  as  a  function  of  both  the  position  and 
speed  of  the  moving  load.  In  their  analysis,  they  approximate 
the  deflections  of  actual  gear  teeth  due  to  moving  loads  by 
modelling  the  tooth  as  a  tapered  Timoshenko  beam.  They  present 
plots  of  normalized  centerline  deflections  for  different  moving 
load  speeds,  and  claim  that  dynamic  stiffness  variations  can  be 
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derived  from  their  results.  However,  as  i.  lustrated  later, 
this  claim  turns  out  to  be  false. 

In  order  to  make  the  theoretical  developments  of  models 
used  to  determine  the  dynamic  response  of  gear  teeth  more  prac¬ 
tical,  some  assumptions  are  made.  One  such  assumption  made  by 
the  first  three  authors  presented,  is  that  the  mass  of  the  gear 
tooth  compared  to  the  gear  body  is  small  and  can  be  neglected. 
Literature  gives  no  hints  to  whether  this  assumption  has  been 
thoroughly  investigated. 

(1.2)  Problem  Statement 

In  this  study,  two  basic  problems  are  investigated.  The 
first  phase  is  to  determine  whether  the  dynamic  response  of  a 
single  spur  gear  tooth  is  dependent  on  the  speed  of  a  moving 
load  acting  on  the  tooth. 

The  second  phase  is  afr^tnv^tigation  to  determine  the 
significance  of  omitting /thel  inertia  of /the  gear  tooth  from  the 
dynamic  deflection  modelV due  so  the  sptall  mass  relative  to  the 
gear  body.  ^ - 

(1.3)  Scope  of  Work 

A  model  based  on  involute  geometry  is  developed  to  automa¬ 
tically  generate  a  spur  gear  tooth  profile  and  finite  element 
mesh,  including  the  rim,  using  a  minimum  of  input  parameters. 
This  model  is  then  used  to  determine  the  effects  of  the  speed 
of  a  moving  load  on  the  deflection  of  a  single  gear  tooth.  Two 
constraint  configurations  are  tested;  one  where  only  the  invo¬ 
lute  profile  and  fillet  regions  are  allowed  to  deform,  the 
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other  with  the  entire  rim  included.  The  results  are  first 
represented  as  normalized  deflections  of  the  tooth  centerline. 
Then  the  tooth  tip  deflection  time  histories  are  studied  for 
the  entire  load  cycle. 

The  second  phase  of  the  work  is  to  model  a  meshing  gear 
tooth  pair  using  two  cantilever  beams  attached  to  moveable 
foundation  masses.  Relative  displacements  of  the  foundation 
masses  as  well  as  beam  deflections  are  determined  for  moving 
load  speeds  bracketing  the  system  resonant  frequency. 


2 .  MODEL  DEVELOPMENT 


In  order  to  effectively  perform  a  static  and  dynamic  ana¬ 
lysis  of  spur  gear  teeth  using  finite  element  techniques,  a 
model  is  needed  to  automatically  generate  a  tooth  profile  and 
the  accompanying  finite  element  mesh  for  different  size  gear 
teeth.  Also,  the  geometry  of  the  tooth  should  be  defined  using 
a  minimum  of  parameters  corresponding  to  those  most  generally 
specified  when  generating  a  tooth  profile.  One  such  list  of 
parameters  is: 


Pressure  Angle  =  0p 

Pitch  Radius  =  RP 

Addendum  =  AD 

Dedendum  =  DED 

Circular  Pitch  =  CIRP 

Backlash  =  BACKL 

Fillet  Radius  =  RF 

Rim  Thickness  =*  RTH 


With  these  parameters,  the  profile  of  any  spur  gear  tooth  can 
be  generated  including  the  rim. 

In  the  proceeding  sections,  the  equations  necessary  to 
construct  the  tooth  profile  using  the  preceeding  parameters  are 
developed,  including  the  implementation  of  these  relationships 
in  a  profile  generation  algorithm.  The  topic  of  finite  element 
mesh  generation  is  also  discussed,  along  with  an  overview  of 
the  mesh  generation  algorithm  used  to  generate  the  grid. 


Later  in  this  chapter,  a  brief  discussion  of  the  plane 
strain  finite  element  type  used  to  model  the  gear  is  included. 
Also,  a  general  treatment  of  a  linear  quadrilateral  element  is 
used  to  help  develop  equations  describing  the  moving  loads  used 
on  the  gear  teeth.  These  relations  are  then  implemented  in  a 
moving  load  generation  algorithm  using  idealized  load  time 
history  equations  for  a  spur  gear  tooth. 

(2.1)  Profile  Generation 

The  profile  generation  sequence  is  divided  into  three 
sections;  determining  relationships,  first  for  the  involute, 
then  for  the  fillet,  and  finally  for  the  rim. 

(2.1.1)  Involute  Generation 

An  involute  curve  is  generated  by  unwrapping  an  inexten- 
sible  cord  from  a  cylinder.  Figure  2-1  illustrates  that  as  the 
cord  is  unwrapped  from  the  cylinder,  point  B  on  the  cord  traces 
an  involute  curve  AC.  The  radius  of  curvature  of  the  involute 
varies  continuously,  being  a  zero  at  point  A  and  increasing 
towards  C.  At  the  instant  shown,  the  radius  is  equal  to  BE,  as 
point  E  corresponds  to  the  instantaneous  center  of  rotation 
about  point  B.  When  generating  the  involute  of  a  spur  gear 
tooth,  the  cylinder  from  which  the  cord  is  unwrapped  corresponds 
to  the  base  circle.  This  concept  is  further  developed  as  shown 
in  Figure  2-2.  Here  the  local  coordinate  system  X-Y,  fixed  to 
the  hub  at  the  base  circle,  is  used  to  determine  relative  coor¬ 
dinates  along  the  involute.  The  parameters  shown  are;  the  base 
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circle  radius  RB,  the  pitch  circle  radius  RP,  the  roll  angle 
and  the  pressure  angle  0p.  The  pressure  angle  is  defined 
by  drawing  a  line  perpendicular  to  the  base  circle  and  passing 
through  the  point  P.  This  line  corresponds  to  the  pressure 
line  or  line  of  action  of  forces  between  meshing  gear  teeth. 
The  point  P  being  the  pitch  point.  From  the  triangle  OBPO,  the 
base  circle  radius  is  defined  in  terms  of  the  pitch  radius  as: 

RB  *  RPcosep  (2-1) 

In  order  to  generate  discrete  points  along  the  involute, 
8i  is  used  in  place  of  0^  and  allowed  to  vary  from  zero  to  a 
maximum,  0max  *^9i?  i*l,2,3, . . . ,n,  corresponding  to  the  desired 
height  of  the  involute,  as  shown  in  Figure  2-3.  The  maximum 
value  of  6  ,  corresponding  to  the  point  Bn  on  the  tip  of  the 
tooth,  is  found  by  writing  the  equation  for  the  triangle 
OAnBnO; 

<RB0max>2  +  Rb2  =  RO2 
Solving  for  0max  gives; 


0. 


max 


=  ^RQ2-RB2  j  ** 


(rad) 


(2-2) 


where  RO  is  the  outer  radius  defined  as  the  sum  of  the  pitch 
radius  and  the  addendum.  Each  increment  of  9 i  produces  a  point 
on  the  involute  progressively  further  from  the  base  circle.  In 
terms  of  the  local  axis  system,  X-Y,  the  coordinates  of  the 
points  Bi  are  determined  from  the  geometry  shown  in  Figure 


Figure  2-3:  Construction  Of  an  involute  curve 


2-4.  In  simplified  form  the  equations  for  the  X  and  Y  coor¬ 
dinates  are: 


XBi  *  -RB(sin8i  "  QicosSi)  (2-3) 

YBi  *  RB ( cos  6i  +  SisinQi  -  1)  (2-4) 


where  RB0i  is  the  arc  length  from  the  origin  of  the  X-Y  system 
to  point  Ax  * 

Next,  those  equations  defining  the  overall  geometry  or 
size  of  the  tooth  are  presented.  From  Figure  2-5,  it  can  be 
seen  that; 

(RB0P  )2  +  RB2  -  RP2 


from  which; 


(rad) 


Also  from  Figure  2-5,  it  is  obvious  that; 


which  can  also  be  written  as; 


(2-5) 


ep  -  tan-tyP)  (2-6) 

where  ©P  is  expressed  in  radians.  During  the  process  of 
calculating  actual  tooth  dimensions,  equation  (2-6)  serves  as  a 
useful  derivational  check  on  qP  .  With  ©P  an<*  9P  angle 
4>  is  written  as  the  difference  of  the  two  previous  angles; 

-ep  (2-7) 
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Figure  2-5:  Spur  gear  tooth  geometry 


And  finally#  from  Figure  2-5#  a  is  found  to  be; 

.  -  eP  -  t.n-i(ef )  ♦  £§§*  (j-8) 

where  CIRTH  is  the  circular  thickness  measured  on  the  pitch 
circle#  given  by; 

CIRTH  «  CIR?  -  BACKL  (2-9) 

2 

Since  one  leg  of  <*  passes  through  the  tooth  center#  this  angle 
is  well  suited  for  transforming  the  involute  coordinates  from 
the  X-Y  axis  system  to  a  system  whose  Y  axis  passes  through  the 
center  of  the  tooth#  such  as  Y*  shown  in  Figure  2-6. 

When  analyzing  a  gear  tooth  to  determine  stresses#  deflec¬ 
tions#  etc.#  it  is  very  advantageous  to  make  full  use  of  the 
axisymmetric  properties  of  the  tooth.  The  involute  points 
generated  relative  to  the  X-Y  axis  system  are#  therefore# 
transformed  into  another  system  X'-Y'#  taking  full  advantage  of 
these  properties. 

Using  the  pitch  point  on  one  side  of  the  tooth  as  a 
reference#  as  seen  in  Figure  2-6#  a  vector  r'  is  drawn  from  the 
pitch  point#  Bp#  to  the  gear  center#  0#  which  defines  the  ori¬ 
gin  of  the  X'-Y'  coordinate  system.  The  vector#  r',  is  com¬ 
posed  of  two  vectors#  R  and  r; 

r'  «  R  +  r  (2-10) 


where; 


Figure  2-6:  Definition  of  axisymmetric 
coordinate  system 
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R  ■  dj*  +  RBsinai ' 


(2-11) 


with; 

d  3  RBcosa 

and; 

r  =  xl  +  y J  (2-12) 

(x  and  y  are  the  coordinates  of  point  Bp  calculated  in  terms  of 
X-Y  using  equations  (2-3)  and  (2-4)). 

In  the  new  coordinate  system,  the  coordinates  of  point  Bp 
are  now  defined  as; 

X'  =  RBsinot  +  xcos  a  +  ysina  (2-13) 

Y' .  »  d  -  xsinct  +  ycosa  (2-14) 

Figure  2-7  illustrates  more  clearly  the  elements  comprising 
equations  (2-13)  and  (2-14). 

In  the  profile  generation  algorithm,  included  in  Appendix 
1,  eleven  points  are  calculated  along  the  involute.  Equations 
(2-2),  (2-3),  (2-4),  (2-5),  (2-8),  (2-13),  and  (2-14)  are  used 
directly  to  calculate  the  point  coordinates  in  the  X'-Y'  axis 
system. 

(2.1.2)  Fillet  Generation 

In  the  present  work,  two  different  spur  gear  tooth 
geometries  are  considered;  low  contact  ratio  gearing  ( LCRG ) 
and  high  contact  ratio  gearing  ( HCRG ) .  By  definition,  the  con¬ 
tact  ratio  is  the  length  of  the  path  of  contact  of  mating  gears 
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divided  by  the  base  pitch.  More  practically,  it  can  be  thought 
of  as  the  average  number  of  teeth  in  contact  during  the  meshing 
cycle.  A  high  contact  ratio  gear  is  one  which  has  at  least  two 
teeth  in  contact  at  all  times. 

One  of  the  main  differences  between  the  two  forms  (LCRG  and 
HCRG)  is  the  fillet  transition  (see  Figure  2-8).  For  an  actual 
low  contact  ratio  gear,  the  fillet  radius  is  placed  tangent  to 
the  involute  and  the  root  circle  as  shown  in  Figure  2-8a.  The 
amount  of  overlap  of  the  involute  may  be  different  for  any 
given  tooth  design,  in  the  current  model,  however,  the  fillet 
radius  is  calculated  to  fit  tangent  to  the  involute  at  the  base 
circle  and  tangent  to  the  root  circle,  as  shown  in  the  modified 
case  (see  Figure  2-8b) . 

When  designing  the  high  contact  ratio  gears,  the  fillet 
region  is  undercut  to  provide  additional  clearance  for  the 
engaging  teeth.  Also,  the  HCRG  tooth  is  generally  longer  due 
to  addendum  or  other  profile  modifications,  thus  the  radial 
distance  between  the  base  and  root  circles  is  also  extended  as 
shown  in  Figure  (2-8c). 

Given  the  gear  parameters  defined  for  a  particular  gear, 
the  following  equation  can  be  used  to  determine  whether  the 
gear  is  a  low  or  high  contact  ratio  gear  [10]. 


RR2  +  2RF  RR  >  RB2 


(2-15) 


In  equation  (2-15)  RF  is  the  fillet  radius  specified  for  a 
given  tooth.  If  this  inequality  is  satisfied,  the  tooth  is 
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classified  as  LCRG.  This  means  that  the  specified  fillet  radius 
will  overlap  the  involute ,  and  thus  must  be  changed  to  fit  the 
modified  form  as  described  in  Figure  2-8b. 
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Figure  2-9a  illustrates  the  geometry  used  in  developing 
the  LCRG  fillet  radius  relations.  From  Figure  2-9a  it  is 
obvious  that  the  fillet  radius  needed  to  make  the  transition 
from  the  base  circle  to  the  root  circle  will  have  to  be  larger 
than  the  radial  distance  DD  as  shown.  The  equation  of  the 
given  triangle  is; 


(RF  +  RR)2  -  (RR  +  DD)2  +  RF2 


(2-16) 


By  expanding  and  rearranging  terms,  equation  (2-16)  can  be 


written  for  RF  as; 


2RR  DD  +  DD" 


(2-17) 


This  then  gives  the  equation  of  the  fillet  radius  which  will 
fit  tangent  to  the  involute  at  the  base  circle,  and  tangent  to 
the  root  circle. 

Rewriting  equation  (2-15)  with  the  inequality  reversed 
gives  the  equation  defining  a  HCRG.  For  HCRG  the  transition  is 


RR2  +  2  RF  RR  <  RB2 


(2-18) 


divided  between  a  radial  line  tangent  to  the  involute  and  a 
fillet  radius  from  the  end  of  the  tangent  line  to  the  root 
circle  (see  Figure  2-8c).  Instead  of  calculating  a  new  fillet 
radius,  as  done  for  LCRG,  equation  (2-16)  is  used,  along  with 
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the  specified  fillet  radius ,  to  calculate  the  length  of  the 


radial  line  DR  (see  Figure  2-9b).  Rewritten  in  another  form, 
equation  (2-16)  becomes; 

DD2  +  2RR  DD  -  2RF  RR  «  0  (2-19) 

and  can  be  used  to  determine  the  radial  distance  DD  spanned  by 
the  fillet  radius.  Using  the  positive  root  of  the  quadratic 
equation  (2-19)  for  DD  yields; 

DD=  -2RR+  (4RR2  ♦  8RF  RR) : *  ,2_20) 

DD  is  then  subtracted  from  the  difference  between  the  base  and 
root  radii  to  give  the  length  of  the  radial  tangent  line. 

DR  «  (RB  -  RR)  -  DD  (2-21) 

When  programming  the  preceeding  equations  to  calculate  the 
fillet  coordinates,  eight  equally  spaced  points  are  used.  For 
LCRG,  the  arc  AOB  is  divided  up  into  eight  equal  angles,  6j[ 
(see  Figure  2-10a).  Coordinates  of  successive  points  are 
calculated  by  adding  di's  together  for  i>l,2...,8  until  the 
arc  from  A  to  B  is  generated.  The  coordinates  of  Bj.  in  Figure 
2-10 b  are  found  from; 


xBi  ■  xa  "  RFcose 
Yfli  -  YA  -  RFsine 


For  HCRG,  the  radial  distance  required  for  the  fillet 
radius,  and  the  radial  distance  of  the  tangent  line  may  vary 


from  one  gear  design  to  another.  To  insure  equal  point 
spacing,  integer  arithmetic  is  used  to  weight  the  number  of 
points  between  the  radial  portion  and  the  fillet  radius 
according  to  their  respective  sizes.  This  is  done  to  facili¬ 
tate  the  finite  element  mesh  generation  routine  discussed 
later.  As  in  the  involute  profile  generation  scheme,  the 
points  for  the  fillet  are  calculated  in  the  X-Y  coordinate 
system,  and  then  transformed  to  the  X'-Y'  system. 

(2.1.3)  Rim  'Generation 

With  the  involute  and  fillet  defined,  the  rim  is  then 
generated.  As  stated  previously,  the  fillet  radius  is  placed 
such  that  it  is  tangent  to  the  root  circle  for  both  LCRG  and 
HCRG.  From  this  tangent  point,  the  rim  of  the  gear  is  added  by 
drawing  an  arc  on  the  root  circle.  The  distance  the  arc  is 
extended  on  either  side  of  the  gear  is  approximately  equal  to 
the  circular  thickness  of  the  gear  as  shown  in  Figure  2-11. 

The  angle  ADA  is  determined  from; 

ADA  *  sin-1  (XP^CIRTH,  _  sin-1 (||)  (2-22) 

This  angle  is  then  divided  into  six  equal  segments  and  the 
coordinates  of  points  on  the  rim  are  calculated  using  a  method 
similar  to  that  shown  in  Figure  2-10.  From  the  last  point  on 
the  root  circle,  coordinates  for  a  radial  line  extending  inward 
a  distance  equal  to  the  specified  rim  thickness  RTH  are  calcu¬ 
lated.  To  complete  the  tooth  profile,  coordinate  points  on  the 
inner  portion  of  the  rim  are  calculated  using  a  technique  simi¬ 
lar  to  that  used  for  the  outer  rim. 
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(2.2)  Finite  Element  Mesh  Generation 

No  absolutely  correct  method  has  been  found  to  model  a 
system  using  a  finite  element  mesh,  even  though  the  topic  of 
mesh  development  has  been  treated  quite  extensively.  With  dif¬ 
ferent  element  types  and  solution  techniques,  several  equally 
valid  methods  are  available  for  any  particular  application. 

Indeed,  Cook  [14]  states  that  although  an  optimum  mesh  can 
be  determined  by  requiring  that  element  boundaries  follow  lines 
of  constant  strain,  this  optimum  condition  only  exists  for 
one  set  of  loading  conditions.  As  the  load  changes,  so  does 
the  optimum  mesh  configuration,  and  for  problems  involving 
other  than  static  loading,  the  difficulties  are  compounded. 

However,  when  developing  a  mesh,  simple  guidelines  can  be 
followed  which  will  produce  a  well  enough  refined  mesh  to 
obtain  more  than  satisfactory  results.  To  mention  a  few;  ele¬ 
ment  boundaries  should  be  aligned  with  structural  or  geometric 
boundaries  and  principal  load  trajectories,  element  aspect 
ratios  should  be  kept  low  (less  than  7),  and  when  different 
element  sizes  are  used  transitions  between  different  size  ele¬ 
ments  must  be  gradual  (mesh  grading). 

The  finite  element  mesh  generation  algorithm  used  for  this 
analysis  was  developed  in  accordance  with  the  preceeding  rules, 
as  well  as  maintaining  computational  efficiency.  Figures  2-12 
and  2-13  show  the  nodes  and  elements,  respectively,  for  a  low 
contact  ratio  gear,  and  Figures  2-14  and  2-15  illustrate  the 
high  contact  ratio  geometry  and  also  the  varying  rim  thickness. 

The  grid  consists  of  319  nodes  and  276  quadrilateral  ele- 
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Figure  2-13:  Elements 
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Figure  2-14:  (HCRG) 
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Figure  2-15:  (HCRG) 
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ments.  Ten  equally  divided  vertical  rows  are  used  to  form  the 
involute  portion  of  the  gear  (elements  1-100).  Nodes  on  the 
surface  of  the  involute  correspond  to  the  actual  coordinate 
points  calculated  in  the  profile  generation  routine.  Close  to 
the  surface  of  the  involute,  the  element  spacing  is  small  pro¬ 
viding  additional  stiffness  for  the  application  of  the  load. 
Towards  the  center  of  the  tooth  the  element  spacing  is  greater 
where  less  stiffness  is  needed. 

The  transition  from  the  end  of  the  involute  to  the  root 
circle  is  accomplished  using  one  of  the  two  techniques 
described  in  section  (2.1.2).  For  both  LCRG  and  HCRG,  eight 
equally  spaced  rows  of  elements  are  used  for  the  transition, 
again  using  the  actual  coordinate  points  calculated  in  the  pro¬ 
file  generation  section  as  surface  nodes.  When  using  the  HCRG 
transition,  with  the  radial  line  and  fillet  radius,  the  eight 
surface  nodes  are  divided  between  the  two  sections  keeping 
nodal  spacing  as  even  as  possible. 

In  order  to  maintain  continuity  between  different  gear 
geometry  finite  element  meshes,  elements  1  through  204  remain 
the  same  size  relative  to  the  actual  tooth  sizes.  In  other 
words,  no  changes  are  made  in  the  grid  geometry  during  the 
generation  of  a  particular  gear  model.  The  exception  to  this 
rule  is  that  elements  205  through  276  do  vary  in  size  depending 
on  the  rim  thickness.  Figures  2-14  and  2-15  show  this  variation. 

The  algorithm  containing  the  equations  developed  for  the 
profile  geometry,  as  well  as  those  relationships  used  to  create 
the  finite  element  mesh  is  included  in  Appendix  1. 


(2.3)  Element  Description 
(2.3.1)  Planar  Elements 


Two  element  types  are  considered  for  this  analysis;  plane 
strain  and  plane  stress.  Due  to  the  geometry  and  loading  con¬ 
ditions  of  the  tooth,  it  is  modelled  as  a  plane  elastic  problem. 
A  plane  body  is  a  region  of  uniform  thickness  contained  within 
two  parallel  planes.  When  the  thickness  of  the  body  is  large 
compared  to  the  lateral  dimensions,  the  problem  is  considered  to 
be  plane  strain.  If  the  thickness  is  small,  it  is  considered  to 
be  plane  stress.  The  difference  between  plane  strain  and  plane 
stress  elements  is  evidenced  in  the  material  property  matrices. 
For  isotropic  materials,  the  material  property  matrix  for  the 
case  of  plane  strain  is; 
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where  E«30.E6  is  the  elastic  modulus  and  y*0.3  is  Poisson's 
ratio.  The  matrix  multiplication  factor  is  larger  for  plane 
strain  than  for  plane  stress. 


£(I-m) 

plane  strain:  (|-+m)(I-2m)  “  4.0385E7 
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plane  stress: 
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=  3.2967E7 


The  matrix  elements  are  also  larger  (except  for  element  3,3)  for 
plane  strain.  When  combined  with  the  strain  displacement  rela¬ 
tions  to  form  the  stiffness  matrix,  these  differences  result  in 
an  increase  in  the  stiffness  for  plane  strain  compared  to  plane 
stress . 

The  thickness  of  the  tooth  used  in  the  analysis  is  0.25 
inches.  Comparatively,  the  largest  and  smallest  planar  dimen¬ 
sions  on  the  actual  tooth  (not  including  the  rim)  are  0.224  and 
0.081  inches,  respectively.  Based  on  the  dimensions  it  is  dif¬ 
ficult  to  make  a  judgement  on  the  correct  element  type  for  this 
analysis . 

Figure  2-16  shows  representative  static  deflection  curves 
for  the  plane  strain  and  plane  stress  element  types.  The  addi¬ 
tional  stiffness  of  the  plane  strain  element  is  noted.  Since  the 
difference  in  deflections  between  the  two  element  types  is  small, 
the  plane  strain  element  type  is  chosen  for  this  analysis. 

(2.3.2)  General  Element  Description 

The  plane  strain  element  described  earlier  can  be  repre¬ 
sented  by  a  linear  quadrilateral  element  similar  to  that  shown 
in  Figure  2-17a.  The  intersection  of  the  lines  which  bisect 
the  sides  of  the  element  form  a  normalized  coordinate  system 
Sn  ,  where; 


.2  .4  .6  .8  1.0  .2  .4  .6  .8  1.0 

X/L  X/L 

-  PLANE  STRAIN, - PLANE  STRESS 

STATIC  DEFLECTION  CURVES:  COMPARISON  OF  DIFFERENT  ELEMENT  TYPES 

Figure  2-16: 
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Between  element  corners,  £  and  n  vary  from  -1  to  +1.  The 
displacements  within  the  elements  can  be  written  in  terms  of 
the  shape  functions  for  each  node,  Ni,  Nj,  Nk  and  Nl  as; 

U  *  UiNi (x,y)  +  UjNj (x,y)  +  UkNk(x,y)  +  UlNl(x,y)  (2-23) 

where  Ui,  Uj,  etc.  define  the  magnitudes  of  the  displacements. 
If  all  nodal  displacements  are  zero  except  for  the  coefficient 
of  Ni(x,y),  which  is  defined  as  unity,  the  displacement  from 
node  i  to  the  other  nodes  will  decrease  from  unity  to  zero. 
Using  the  parameters  shown  in  Figure  2-17b  the  shape  function 
for  node  i  going  from  i  to  j  is; 


Ni (x,y)  -  (hzL) 


(2-24) 


where  L  is  the  length  between  nodes  i  and  j  in  the  direction 
of  5. 


(2.3.3)  Moving  Loads 

An  arbitrary  load,  P(  £,  t  ),  normal  to  £  is  introduced 
whose  components  are;  Px(£,t),  Py(  £,t)  (see  Figure  2-17bl  The 
effect  of  the  force  P(£,t)  on  node  i  can  be  represented  by  the 
integral  of  the  load  times  the  shape  function  and  thickness 
in  the  direction  of  from  0  to  L.  Component  wise; 
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where;  At  is  the  element  thickness  (assumed  to  be  unity). 
Inserting  the  shape  function  for  node  i  into  equations  (2-25) 
and  (2-26)  yields; 


h 


Fx1  « 
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J'  PxU,t)  <^)d£ 

(2-27) 

Fyi  * 

f* Py<S,t)  <^)ds 

(2-28) 

0 

Equations  (2-27)  and  (2-28)  give  the  load  history  at  node  i  as 
a  function  of  £(t),  resulting  from  the  arbitrary  load  P(£,t). 

Conversely,  the  load  history  at  node  j  is  determined  by 
considering  the  shape  function  obtained  when  going  from  node 
j  to  i  with  j  at  zero  and  i  at  unity.  Here  the  shape  function 
starts  at  zero  and  increases  to  unity  as; 

Nj(x,y)  -  ^  (2-29) 

Substituting  equation  (2-29)  into  equations  (2-25)  and  (2-26) 
yields  the  force  in  the  x  and  y  directions  experienced  at  node 
j,  resulting  from  P(£,t). 
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Equations  (2-27),  (2-28),  (2-30)  and  (2-31)  can  be  used  to 
represent  a  moving  load  by  introducing  the  Dirac  Delta 
Function.  When  used  in  an  integral,  it  translates  a  given 
function  to  the  origin  and  gives  the  value  of  a  function  at  a 
given  time  at  the  origin.  The  argument  of  the  delta  function 
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takes  the  form  of  the  variation  in  the  position  variable.  For 
a  moving  load  with  constant  velocity,  the  change  in  position 
is  given  by  the  velocity  times  the  time.  The  arbitrary  moving 
load  then  takes  the  form; 

P(£,t)  *  P(t)«s<5-Vot> 

where  Vo  is  the  velocity  and  t  the  time.  Using  the  delta 
function  in  the  integrand  results  in  all  occurences  of  being 
replaced  by  Vot.  Thus,  the  four  force  equations  become; 


Fx1 

»  Px(t) 

t-* 

i 

(2-32) 
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*  Py( t) 
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L 

(2-33) 

FxD 

*  Px(t) 

<^> 

(2-34) 

Fy3 

*  Py(t) 

1 

(2-35) 

Plotting  these  equations  as  a  function  of  time  where  the  magni¬ 
tude  of  P(t)  is  constant,  yields  to  general  force  histories 
(see  Figure  2-18)  for  a  load  moving  from  i  to  j. 


Figure  2-18:  Linear  force  histories 
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For  an  actual  meshing  gear  set,  the  speed  of  the  moving 
load  on  a  single  tooth  is  not  constant,  but  varies  linearly 
with  time.  The  time  varying  speed  can  be  seen  to  be  (see 
Appendix  2) ; 

V(t)  ■  RB  m2 t 


Now,  the  force  equations  take  a  different  form  with  3  being 
replaced  by  the  displacement  resulting  from  the  above  velocity; 


S(t) 


RBu2t 

2 


where  A  can  replace  the  quantity  RB*2/2; 


S(t)  -  At2 


For  the  time  varying  load  the  force  equations  then  become; 


F*  -  P(t) 

<^2» 

(2-36) 

Fj  -  P(t) 

(2-37) 

where  the  x  and  y  subscripts  are  assumed 


3.  DISCUSSION  OF  NAGAYA  ANALYSIS 


Although  the  problem  of  theoretically  analyzing  dynamic 
gear  tooth  deflections  has  been  treated  extensively  [1-51(10], 
models  addressing  the  problem  assume  that  the  variation  in 
tooth  stiffness  can  be  approximated  using  a  static  deflection 
analysis.  These  models  assume  that  the  gear  hubs  act  as  rigid 
bodies  and  that  the  teeth  act  as  variable  stiffness  springs. 
The  stiffness  of  the  teeth  varies  with  the  contact  position 
along  the  tooth  and  is  generally  arrived  at  using  a  static 
deflection  analysis  such  as  the  one  developed  by  Weber  [12]. 
Recently,  K.  Nagaya  and  S.  Uematsu  [7]  proposed  that  since  the 
contact  point  moves  along  the  tooth  during  the  meshing  cycle, 
the  dynamic  load  response  should  be  considered  as  a  function  of 
both  the  position  and  the  speed  of  the  moving  load.  In  their 
paper  they  generate  plots  of  normalized  gear  tooth  centerline 
deflection  curves  from  which  they  claim  the  equivalent  spring 
constant  of  gear  teeth  can  be  determined. 

(3.1)  Approximating  A  Gear  Tooth  with  a  Timoshenko  Beam 

In  Nagaya' s  analysis,  the  differential  equations  for  a 
taper id  Timoshenko  beam  are  written  and  solved,  in  the  form  of 
an  eigenvalue  problem,  from  which  a  modal  response  analysis  is 
used  to  determine  tooth  deflections  due  to  moving  loads. 
Nagaya  assumed  a  load  of  constant  magnitude  moving  along  the 
beam  at  a  constant  velocity  from  the  tip  to  the  base  of  the 
tooth  (see  Figure  3-1). 

Using  Kara's  [8]  assumption  for  the  profile  of  gear 
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L  -  Ll=  2,25m 
Ho  =  2.48m 
Ll/L  =  ,34 

S  =  (AoLa/Io)*=  4,76 


Figure  3-1:  Tapered  Timoshenko  beam 


teeth,  Nagaya  claimed  that  the  deflections  obtained  using  the 
beam  approximation  were  applicable  to  any  spur  gear  defined  by 
the  parameters 
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where  m  is  the  module,  L,  LI,  Ho  are  shown  in  Figure  3-1,  and  S 
is  the  slenderness  ratio.  The  module,  m,  is  the  pitch  diameter 
divided  by  the  number  of  teeth,  measured  in  inches.  When  ana¬ 
lyzing  a  gear  tooth  the  above  parameters  are  used  to  describe 
the  Timoshenko  beam  used  for  the  approximation.  An  example  of 
such  a  comparison  is  shown  in  Figure  3-2  where  the  approxi¬ 
mating  Timoshenko  beam  is  shown  superimposed  onto  the  gear 
tooth  used  in  the  finite  element  analysis  of  Chapter  4. 
Instead  of  the  beam  lying  tangent  to  the  involute  of  the  actual 
test  gear  as  shown,  it  should  have  passed  through  the  tip  of 
the  involute  as  illustrated  by  the  inset  figure.  The  inset  is 
a  correct  representation  of  Kara's  assumption  for  the  profile 
of  spur  gear  teeth.  This  discrepancy,  is  solely  attributable 
to  the  use  of  backlash  when  defining  the  gear  geometry. 
Backlash  effectively  decreases  the  width  of  the  tooth.  In 
order  to  better  compare  the  finite  element  analysis  to  Nagaya 's 
work,  the  tootu,  when  analyzed,  is  constrained  so  that  only  the 


gear  teeth 


portion  void  of  interior  elements  is  allowed  to  deform  (see 
Figure  3-2).  The  foundation  and  rim  are  constrained  against 


motion.  Later  on,  the  deflection  of  the  gear  tooth  is  again 
analyzed  with  the  tooth,  foundation,  and  rim  allowed  to  deform. 

(3.2)  Interpretation  of  Nagaya  Results 

When  presenting  his  findings,  Nagaya  plotted  normalized 
tooth  centerline  deflections  versus  normalized  load  position 
for  different  moving  load  speeds.  Figure  3-3,  taken  directly 
from  reference  [7],  illustrates  these  results.  The  solid  cur¬ 
ves  in  the  figure  represent  normalized  tooth  centerline  deflec¬ 
tions,  each  one  for  a  different  normalized  velocity,  V*.  The 
vertical  arrows  labelled  T  represent  the  load  position  relative 
to  the  length  of  the  tooth,  where  X/L  is  the  ratio  of  the  posi¬ 
tion  of  the  load  on  the  tooth  relative  to  the  total  length  L. 
The  dotted  lines  are  the  static  curves  obtained  from  the  Karas 
analysis.  By  normalizing  these  parameters;  deflection,  load 
position,  and  velocity,  the  results  then  become  applicable  to 
any  size  gear  tooth. 

Non-dimensional  deflections  can  be  represented  by; 

W*  »  AoEW/PL 

where 

Ao  *  Area  of  the  base  of  the  tooth  (in^) 

E  *  Elastic  modulus  (Psi) 

W  *  Actual  tooth  deflection  in 

direction  of  applied  load  (in) 

P  »  Applied  load  (lbs) 

L  *  Extended  tooth  length  (in) 
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Dynamic  deflection  curvet  for  the  gear  teeth  with  various 
moving  loads 


Figure  3-3 


The  velocity  in  normalized  form  is  written  as; 


V*  -  V/  JWp 

where; 

V  3  Speed  of  moving  load  (in/sec) 
E  3  Elastic  modulus  (Psi) 
p  *  Material  density  (lb/in^) 


In  this  equation  the  denominator  represents  the  wave  velocity 
in  bars.  Finally,  the  position  of  the  moving  load  is  given  by; 


T  «  Vt/ (L-Ll ) 

where; 

V  3  Speed  of  the  moving  load  (in/sec) 
t  *  Elapsed  time  (sec) 

( L-Ll )  *  Actual  tooth  height  (in) 


From  the  plots  shown  in  Figure  3-3,  Nagaya  claims  that  the 
deflections  of  gearteeth,  subjected  to  moving  loads,  vary  with 
the  speed  of  the  moving  load.  That  is,  for  the  same  values  of 
T,  the  displacements  are  directly  related  to  the  speed  of  the 
load.  He  states  that  for  slowly  moving  loads,  the  dynamic 
response  reduces  to  the  case  of  a  step  function  impact  load  for 
small  values  of  T  (see  T»0.1,  V**0.001  in  Figure  3-3).  Since 
Figure  3-3  indicates  that  the  dynamic  response  is  dependent  on 
the  moving  load  speed  (due  to  effects  of  inertia  forces  of  the 
mass  of  the  tooth),  Nagaya  states  that  the  stiffness  of  the 
tooth  must  also  depend  on  the  moving  load  speed.  he  then 


claims  that  the  varying  tooth  stiffnesses  can  be  determined 


from  these  plots.  However,  as  later  demonstrated,  Nagaya’s 
claim  that  the  response,  and  therefore  the  stiffness,  is  depen¬ 
dent  on  the  speed  of  the  moving  load  is  a  false  one. 

A  major  portion  of  the  present  work  is  directed  towards 
substantiation  of  this  conclusion. 


4.  FINITE  ELEMENT  ANALYSIS 


The  deflections  of  single  spur  gear  teeth  with  moving 
loads  acting  on  them  are  determined  using  finite  element  analy¬ 
sis.  A  single  gear  tooth  is  used  for  six  different  moving  load 
cases.  First,  the  same  moving  load  scheme  used  by  Nagaya  [71 
(constant  magnitude  and  speed)  is  applied  to  the  tooth,  which 
is  constrained  according  to  the  Timoshenko  beam  approximation. 
Then  the  load'  application  on  the  tip  of  the  tooth  is  changed 
slightly  and  the  test  repeated  on  the  tooth  with  the  same 
constraints.  The  two  preceeding  load  cases  are  then  applied  to 
a  tooth  allowing  the  entire  model  to  deform,  including  the  rim. 
Finally,  an  idealized  load  function,  with  variable  load  magni¬ 
tude  and  speed,  is  applied  to  the  tooth  using  both  constraint 
cases. 

(4.1)  Description  of  Test  Gear 

The  gear  used  as  the  model  for  this  analysis  was  selected 
at  random.  The  parameters  used  to  define  the  geometry  of  the 
gear  are? 


ep 

- 

Pressure  angle 

20® 

RP 

* 

Pitch  radius 

- 

1.75 

(in) 

AD 

* 

Addendum 

- 

0.125 

(in) 

DED 

* 

Dedendum 

* 

0.175 

( in) 

CIRP 

* 

Circular  pitch 

- 

0.3927  (in) 

BACKL 

* 

Backlash 

- 

0.01 

(in) 

RF 

a* 

Fillet  radius 

. 

0.05 

(in) 

RTH  =  Rim  thickness  *  d.6  (in) 

At  -  Tooth  thickness  ■  0.25  (in) 

Figure  4-1  shows  the  finite  element  model  of  the  tooth. 
The  test  gear  is  a  low  contact  ratio  gear  (contact  ratio  * 
1.74) . 

(4.2)  Determination  of  Normalized  Plotting  Parameters 
and  Their  Application  to  the  Gear  Tooth 

As  stated  previously,  the  normalized  deflections  of  the 
gear  tooth  are  calculated  using  those  parameters  specified  in 
Kara's  assumption  for  the  profile  of  gear  teeth.  Thus,  when 
the  deflections  are  plotted,  the  only  term  in  the  normalized 
deflection  equation  taken  directly  from  the  gear  analysis  is 
the  deflection  of  the  tooth  centerline  in  the  direction  of  the 
applied  load  which  is  perpendicular  to  the  centerline  of  the 
tooth . 

In  Chapter  3.1  the  equations  needed  to  define  the  tooth 
profile  approximation,  according  to  Karas,  are  given.  The  phy¬ 
sical  dimensions,  length,  area,  etc.  are  defined  in  terms  of 
the  module,  m.  For  a  standard  spur  gear  the  module  is  defined 
as  the  pitch  diameter  per  tooth  measured  in  inches,  and  is 
usually  represented  by  the  inverse  of  the  diametral  pitch; 

Module  -  M  -  1/DP  (in)  (4-1) 

where  the  diametral  pitch  is; 

diametral  pitch  ■  DP  ■  n  ^  M  *  8  (4-2) 

CIRP  .  3927 
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Figure  4-1:  Test  gear 


To  further  define  the  test  gear,  the  number  of  teeth  can  be 
calculated  from; 


number  of  teeth  ■  N  *  2RP*DP  =  2 (1.75) (8)  *  28  (4-3) 


Given  either  the  diametral  pitch  or  the  number  of  teeth,  the 
module  can  easily  be  obtained. 


m  ■  0.125  ( in) 


Using  the  value  for  the  module  and  the  relations  of  Chapter 
3.1,  the  dimensions  of  the  approximating  Timoshenko  beam  are 


determined.  The  height  of  the  corresponding  beam  becomes; 


( L-Ll )  *  2.25m  *  0.28125  (in) 


(4-4) 


and  the  extended  length; 


L  *  (£~i)  -  0.42614  (in) 

0*66 


(4-5) 


At  the  base,  the  beam  thickness  is; 


Ho  =  2.48m  »  0.31  (in) 


and  thus  the  area  at  the  base; 


Ao  -  HoAt  ■  (0.31X0. 25)  -  0.775  (in) 


(4-6) 


Given  the  above  parameters,  the  non-dimensional  deflections  can 


be  plotted  using; 


W*  -  AoEW/PL 


(4-7) 


Again,  it  should  be  emphasized  that  although  the  beam  approxi¬ 
mation  (shown  in  Figure  3-2)  does  not  match  the  tooth  exactly. 


V-' v'V 


the  dimensions  just  determined  in  equations  (4-5)  and  (4-6)  are 
still  used  in  the  normalized  deflection  equation  when  plotting 
results  for  the  actual  tooth. 


Nagaya  defines  the  normalized  velocity  of  the  moving  load 
to  be; 


V*  -  (4-8) 

and  plots  four  speeds  corresponding  to  V*  equal  to;  0.01* 
0.005,  0.003,  and  0.001.  Using  equation  (4-8),  the  actual 

velocities,  V,  are; 


- V* - 

V( in/ sec) 

0.01 

2025 

0.005 

1012 

0.003 

607 

0.001 

203 

Table  4-1:  Actual  Velocities 

When  the  load  moves  along  the  involute  at  a  constant  speed,  the 
values  shown  in  the  preceeding  table  are  used  directly. 
However,  for  a  meshing  gear  set,  the  velocity  along  the  invo¬ 
lute  changes  from  zero  to  a  maximum  velocity,  Vmax,  according 
to  equation  (A2-6).  When  the  speed  of  the  moving  load  is 
modelled  using  this  relation,  the  maximum  velocity  is  defined 
to  be  the  velocity  given  in  Table  4-1.  Therefore,  the  speed 
starts  at  zero  and  increases  to  a  maximum  speed  corresponding 
to  those  given  in  the  table.  By  choosing  the  previously  deter¬ 
mined  velocities  of  Table  (4-1)  to  occur  at  the  base  circle 
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where  the  velocity  is  maximum,  the  RPM  of  the  gear  for  each 


non-dimensional  velocity  can  be  determined.  From  equations 
(A2-6)  and  (A2-7)  the  load  cycle  time  for  a  meshing  tooth  can 
be  shown  to  be; 


TF  »  28/Vmax  (4-9) 

where  S  is  the  distance  along  the  involute  from  the  tip  to  the 
base  circle,  and  Vm^v  is  the  velocity  at  the  base  circle. 
Equation  (A2-6)  can  then  be  used  to  calculate  the  angular  velo¬ 
city  of  the  gear; 

Vmax 

OMEGA  *  TF*RB  (rad/sec)  (4-10) 

For  each  non-dimensional  velocity,  load  cycle  times,  TF,  and 
the  RPM's  of  the  test  gear  are  found  to  be; 


V* 

TF ( sec ) 

RPM 

0.001 

2025 

. 2436E-3 

21470 

iW  <T  <T.~MitIll 

1012 

. 4875E-3 

10730 

0.003 

607 

. 8127E-3 

6435 

0.001 

203 

. 2436E-2 

21*7 

Table  4-2:  Variable  Velocity  Parameters 


To  reiterate,  the  times  TF  included  in  Table  4-2  are  those  for 
which  the  velocity  starts  at  zero  at  the  tip  and  increases 
linearly  to  a  maximum  value,  Vmax-  For  a  load  moving  with 
constant  velocity,  the  time  for  the  load  to  move  over  the 
involute  is  simply  the  distance,  S,  divided  by  the  velocity. 

Returning  now  to  the  non-dimensional  parameters  plotted  by 


Nagaya,  for  a  constant  speed  moving  load  the  position  of  the 
load  along  the  involute  is  given  by; 


T  “E=EI 

where  T  varies  from  0.0  at  the  tip,  to  1.0  at  the  root  circle. 
A  value  T*0.8  correponds  to  a  point  near  the  base  circle 
radius  between  nodes  110  and  121  (see  Figure  2-12).  Equation 
(4-11)  is  valid  only  for  constant  speeds,  V. 

(4.3)  Description  of  Dynamic  Loading  Cases 

The  dynamic  deflections  of  single  spur  gear  teeth  are 
generated  using  three  loading  cases;  a  constant  speed  constant 
magnitude  load  with  impact  engagement,  a  constant  speed 
constant  magnitude  load  using  a  finite  load  engagement  rise 
time  (for  these  two  loading  cases  the  load  is  applied  normal 
to  the  tooth  centerline),  and  a  load  with  varying  speed  and 
magnitude.  In  this  last  case  the  load  is  applied  normal  to 
the  involute. 

The  first  of  these  three  loading  cases  is  designed  to  imi¬ 
tate  exactly  the  forcing  function  used  by  Nagaya.  At  time 
equal  to  zero,  a  load  of  1000  lbs,  simulating  an  impact  load, 
is  applied  to  the  tip  of  the  tooth,  and  maintained  until  the 
end  of  the  load  cycle  (i.e.  from  the  tip  to  the  base  circle). 
(Whenever  the  terms  "impact  loading"  are  used,  the  author  is 
describing  a  step  function).  To  simulate  this  loading  con¬ 
dition  for  the  finite  element  analysis,  time  functions  repre¬ 
senting  nodal  load  histories  are  calculated  for  each  node  on 


the  involute.  For  the  load  moving  between  nodes  i  and  j,  the 
force  histories  are  described  by  equations  (2-36)  and  (2-37); 


Fi  -  p(t) 

(4-12) 

Fj  -  P(t)  (2£) 

(4-13) 

where;  L  is  the  distance  between  nodes,  V  is  the  velocity  of 
the  moving  load,  and  t  the  time.  With  several  load  value  data 
points  defined  along  the  involute,  the  finite  element  code  uses 
these  points  and  linearly  interpolates  between  them  to  define 
the  time  functions.  The  time  functions  for  this  loading  case 
are  shown  in  Figure  4-2.  In  Figure  4-2  the  node  numbers 
correspond  to  the  first  eleven  nodes  on  the  right  involute  sur¬ 
face  of  the  tooth. 

To  determine  the  effect  of  impact  load  engagement,  another 
test  is  run  using  a  finite  rise  time  for  the  load  on  the  first 
node.  Instead  of  the  load  applied  all  at  once  at  time  equal  to 
zero,  it  starts  at  zero  and  gradually  increases  to  the  maximum 
value  (see  Figure  4-3).  Here  the  magnitude  of  the  load  is  zero 


Figure  4-3;  Finite  rise  time 
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at  t  equal  to  0.0  and  goes  up  to  3/4  of  PmaX  by  PER/8.  The 
rise  time  is  defined  as  a  fixed  fraction  (PER/8)  of  the  time 
function  period.  As  the  speed  of  the  moving  load  increases , 
the  rise  time  decreases.  Therefore,  the  rise  time  for  V**0.001 
is  ten  times  greater  than  for  V*=0.01.  The  time  functions  for 
this  loading  case  are  shown  in  Figure  4-4.  Only  the  first  time 
function  is  different  between  Figures  4-2  and  4-4. 

The  last  loading  case  tested  is  one  in  which  the  speed 
and  the  magnitude  of  the  load  vary  with  time.  Wallace  and 
Seireg  [9]  give  idealized  relationships  for  the  magnitude  of 
the  load  on  a  gear  tooth  as  a  function  of  time  and  the  contact 
ratio.  They  are; 


P(t)  »  %  Pmax<l-cos(^)) 


for:  0.0  < t <  aTF 


P(t)  *  Pmax 


for:  aTF<t^(l-a)TF  (4-14) 


P(t)  *  h  Pmax (1- cos  (-- )  for:  ( 1-a )  TF  <  t < TF 

where;  TF  is  the  load  cycle  time,  t  is  the  time  along  the  invo¬ 
lute,  and  a  is  a  factor  dependent  on  the  contact  ratio.  A 
value  of  0.29  for  a  is  used,  corresponding  to  a  contact  ratio 
of  1.56.  The  force  history  described  by  equations  (4-14)  is 
plotted  in  Figure  4-5.  By  applying  equations  (2-36)  and  (2-37) 
with  the  load  replaced  by  equations  (4-14),  the  time  functions 
generated  for  this  loading  case  are  like  those  shown  in  Figure 
4-6.  Note  that  the  time  function  period  decreases  as  the  load 
moves  down  the  involute  due  to  increasing  speed. 

Equation  (4-12)  and  (4-13),  along  with  equation  (4-14)  are 
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used  in  a  time  function  generation  algorithm  which  is  included 
in  Appendix  1. 

(4.4)  Finite  Element  Test  Results 

The  results  contained  in  the  proceeding  sections  were 
obtained  using  the  SAP6  finite  element  code,  implemented  on  the 
UNIVAC  1100/80  computer  facility  at  Michigan  Technological 
University. 

(4.4.1)  Comparison  of  Static  Results 

As  a  preliminary  check  on  the  accuracy  of  the  finite  ele¬ 
ment  analysis  technique  applied  to  gear  teeth,  static- 
deflections  obtained  using  finite  elements  are  compared  to 
those  calculated  by  Nagaya  using  Kara's  assumption  for  the  pro¬ 
file  of  gear  teeth.  Comparisons  are  made  with  and  without  the 
rim  included  in  the  analysis.  Figure  4-7  shows  the  plots  of 
the  normalized  centerline  deflections  obtained  using  Timoshenko 
beam  constraints.  The  dashed  lines  are  the  static  deflections 
calculated  by  Nagaya.  From  the  figure  it  is  apparent  that  the 
Timoshenko  beam  (used  to  produce  the  dashed  lines)  is  stiff er 
than  the  tooth.  Going  back  to  Figure  3-2,  it  is  seen  that  the 
beam  is  considerably  larger  than  the  tooth,  especially  towards 
the  base.  Thus  one  would  expect  the  beam  to  be  stiffer.  As 
the  load  is  applied  closer  to  the  base  of  the  tooth  the  dif¬ 
ference  between  the  static  deflection  curves  becomes  less 
exaggerated.  For  T»0.8  the  centerline  of  the  tooth  actually 
deflects  less  than  the  beam.  The  reasons  for  this  are  not 
completely  clear.  One  possible  explanation,  however,  is  the 


fact  that  as  the  load  is  applied  closer  to  the  base,  the  amount 
of  local  deformation  around  the  point  of  load  application 
increases  due  to  increased  nodal  spacing.  This  causes  the 
tooth  centerline  to  deform  around  the  local  deformation,  thus 
decreasing  the  overall  deflection  of  the  gear  tooth.  (Appendix 
3  includes  the  actual  tooth  in  the  statically  deformed  con¬ 
dition,  illustrating  the  increase  in  local  deformation).  In 
these  figures,  the  compatibility  of  the  element  is  not 
violated.  The  deformation  scale  factor  causes  element  overlap. 

With  the  rim  included  in  the  analysis,  the  centerline 
deflections  are  considerably  more  severe  (see  Figure  4-8). 
The  curves  obtained  by  Nagaya,  represented  by  the  dashed  line, 
are  exactly  those  pictured  on  Figure  4-7.  The  purpose  of  this 
set  of  plots  (Figure  4-8)  is  to  emphasize  the  added  flexibility 
afforded  by  the  rim  material.  (Appendix  3  also  contains  the 
tooth  in  the  deflected  state  with  the  rim  included) . 

(4.4.2)  Modal  Analysis  -  Determination  of  Mode  Shapes  and 
Natural  Frequencies 

In  Nagaya' s  paper,  the  differential  equation  for  the  non- 
dimensional  deflection,  W* ,  is  derived  and  then  solved  numeri¬ 
cally.  The  solution  to  the  differential  equation  (eigenvalue 
problem)  includes  an  infinite  number  of  natural  frequencies 
(eigenvalues)  and  an  infinite  number  of  mode  shapes 
(eigenvectors).  However,  he  included  only  the  first  three 
eigensolutions  in  the  dynamic  response  analysis.  It  should 
also  be  mentioned  that  the  differential  formulation  is  done  for 
transverse  vibration  so  only  bending  modes  are  included. 


X/L  X/L 

- -  FEM  WITH  RIM,  -  KARAS'  STATIC 

STATIC  DEFLECTION  CURVES:  COMPARISON  OF  FEM  AND  KARAS  TECNIQUtS 

Figure  4-8:  Rim  included 
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For  the  present  analysis  the  mode  shapes  and  natural  fre¬ 
quencies  are  determined  from  the  finite  element  model  by 
solving  the  general  equation; 

(IK]  -  AIM]){D}  *  0  for  X-  u)2 

for  pairs  of  x  and  (D>.  These  results  are  then  used  in  a  modal 
response  analysis  to  determine  the  response  of  the  tooth  to  the 
moving  loads.  For  this  analysis  ten  modes  are  included,  with 
both  transverse  and  axial  vibration.  Appendix  3  includes  the 
first  few  mode  shapes  and  natural  frequencies  of  the  tooth  for 
both  constraint  cases.  Again,  note  the  difference  in  flexibi¬ 
lity  between  the  two  models  (with  and  without  the  rim). 

The  first  three  natural  frequencies  from  Nagaya's  work 
are  compared  with  those  found  for  the  actual  tooth.  Both 
bending  and  axial  modes  are  included  in  the  finite  element 
analysis,  so  the  first  three  bending  modes  from’  this  analysis 
are  used  for  comparison  (modes  1,3,4)  (see  Table  4-3). 


NAGAYA  FEM 


Mode 

FREQ  (rad/sec) 

Mode 

FREQ  (rad/sec) 

1 

1.092E6 

1 

5.718E5 

2 

3.764E6 

3 

1.436E6 

3 

9.488E6 

4 

2.650E6 

Table  4-3  Comparison  of  Modal  Results 
As  expected,  the  natural  frequencies  of  the  beam  approximation 
are  somewhat  higher  than  those  of  the  tooth  (beam  constraints), 

68 


I  *  *  * >: 


*  (  t  i':  I’  i  < 


partly  due  to  the  additional  material  toward  the  base  of  the 
beam. 

(4.4.3)  Dynamic  Deflections:  Timoshenko  Beam  Constraints 
(4. 4. 3.1)  Impact  Loading 

Shown  in  Figures  4-9a  and  4-9 b  are  the  normalized  cen¬ 
terline  deflections  obtained  using  the  impact  engagement 
loading  case  (see  Figure  4-2).  Results  are  obtained  for  load 
positions  of  *00.1,0.2, ... ,0.8.  Each  solid  line  represents 
the  non-dimensional  centerline  deflection; 

W*  «  AoEW/PL 

due  to  the  applied  moving  load.  Remember  also  that  the  deflec¬ 
tions  are  plotted  as  a  function  of  position; 


and  not  as  a  function  of  time.  So  for  00.1  the  solid  lines 
show  the  normalized  centerline  deflections  for  the  different 
speeds  with  the  load  one  tenth  the  distance  between  the  tip  and 
the  root.  Remember  also  that  the  time  for  the  load  to  move 
from  OO  to  00.1  is  different  for  all  velocities,  V*.  The 
dashed  lines  correspond  to  the  static  deflections  obtained  for 
the  tooth  with  the  load  in  the  position  shown. 

Initial  examination  of  these  plots  suggests  that  the 
displacements  are  indeed  dependent  on  the  speed  of  the  moving 
load.  For  slow  speeds  and  low  values  of  T  (V**0.001,  T-O.l), 
the  deflection  is  approximately  twice  the  static  as  Nagaya 
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DYNAMIC  DEFLECTION  CURVES  OF  SPUR  GEAR  TEETH  WITH  MOVING  LOADS 

Figure  4-9a:  Timoshenko  beam  constraints, 
impact  load  eng  igrmont 
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DYNAMIC  DEFLECTION  CURVES  OF  SPUR  GEAR  TEETH  WITH  MOVING  LOADS 

Figure  4-9b:  Timoshenko  beam  constraints, 
impact  load  engagement 
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claimed.  However,  these  plots  do  not  give  an  accurate  descrip¬ 
tion  of  the  dynamic  deflections. 

A  much  more  representative  picture  is  obtained  when  the 
actual  deflection  of  the  tooth  tip  is  examined  over  the  load 
cycle  using  small  sampling  intervals,  (  T*0.001  for  V*=0.01). 
Figure  4-10  shows  the  true  time  history  of  the  tooth  tip  as  the 
load  moves  from  the  tip  to  the  base  of  the  tooth.  Instead  of 
the  tooth  being  in  a  particular  deformed  state  at  load  position 
T,  dependent  completely  on  the  speed,  it  actually  oscillates 
about  a  datum  with  constant  amplitude  and  frequency.  By 
assuming  that  the  tip  oscillates  about  the  static  position,  the 
amplitude  of  oscillation  is  approximately  twice  the  static,  and 
is  initiated  by  the  impact  load  at  the  beginning  of  the  load 
cycle.  The  only  effect  the  speed  of  the  moving  load  has  is  to 
change  the  number  of  oscillations  per  load  cycle.  Recall  that 
the  time  for  the  load  to  move  from  T=0.1  to  T=0.8  is  ten  times 
greater  for  V**0.001  than  for  V*=0.01.  Therefore,  approxima¬ 
tely  ten  times  more  oscillations  occur  for  the  slower  of  the 
two  speeds. 

From  Figures  4-9  and  4-10  we  can  then  conclude  that  the 
deflections  of  the  tooth  do  not  depend  so  much  on  the  speed  of 
the  moving  load,  but  on  the  tooth  position  at  the  instant  the 
deflection  sample  is  taken.  By  roughly  lining  up  the  position 
and  deflection  on  the  tip  deflection  curves,  the  position  of 
the  tip  of  the  tooth  shown  in  Figures  4-9a  and  b  can  easily  be 
duplicated. 

The  datum  about  which  the  tooth  oscillates  is  determined 
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Figure  4-10 


by  repeating  the  previous  analysis  with  the  system  critically 
damped.  The  transient  caused  by  the  impact  load  is  then 
"filtered"  out  with  only  the  steady  state  response  remaining. 
Figure  4-11  shows  the  non-dimensional  centerline  deflections 
for  all  four  moving  load  speeds.  From  the  figure  it  is 
apparent  that  all  centerline  deflections  lay  over  the  static 
curve.  To  verify  this  claim,  the  tooth  tip  deflection 
histories  are  again  plotted.  In  Figure  4-12  the  plots  show  the 
tip  following  the  static  curve. 

(4. 4. 3. 2)  Finite  Engagement  Rise  Time  Loading 

In  this  test  the  tooth  is  subjected  to  the  moving  load 
conditions  illustrated  in  Figure  4-4  where  the  load  on  the 
first  node  is  gradually  applied  over  a  time  of  PER/8.  This 
loading  case  produces  significantly  different  results  compared 
to  the  impact  load  test.  (Since  the  normalized  centerline 
deflection  curves  do  not  accurately  represent  the  dynamic 
deflection  phenomenon,  they  are  not  included).  Figure  4-13 
gives  the  tooth  tip  deflection  history  for  this  loading  con¬ 
dition.  The  tooth  still  oscillates  about  the  static  position 
with  constant  frequency,  but  the  amplitude  varies  significantly 
with  speed.  The  reason  for  this  change,  from  the  previous  load 
case,  is  best  explained  by  again  considering  the  load  engage¬ 
ment  rise  time. 

In  Chapter  4.3,  Figure  4-3,  the  rise  time  is  defined  as  a 
fixed  fraction  of  the  time  function  nodal  period  (PER/8).  For 
V**0.01,  the  rise  time  is  ten  times  less  than  for  V**0.001. 
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It  should  then  be  obvious,  that  as  the  speed  of  the  load 
increases,  the  time  function  for  the  first  node  begins  to 
approximate  the  step  function  impact  load  of  Figure  4-2. 
Examination  of  the  deflection  amplitude  for  V*=0.01  of  Figure 
4-13  shows  it  to  be  nearly  the  same  as  V*=0.01  of  Figure  4-10. 

(4. 4. 3. 3)  Wallace  -  Seireg  Loading 

For  this  loading  case  (time  function  shown  in  Figure  4-6) 
both  the  speed  and  magnitude  of  the  load  vary  with  time.  With 
a  contact  ratio  of  1.56,  the  load  doesn't  reach  the  maximum 
value  of  1000  lbs  until  it  is  between  the  second  and  third 
nodes.  Since  the  magnitude  increases  smoothly  and  gradually, 
no  abrupt  load  changes  are  encountered. 

The  tooth  tip  deflection  history  curves  for  this  loading 
case  are  included  in  Figure  4-14.  Due  to  the  nature  of  the 
speed  variation,  the  deflections  are  plotted  as  a  function  of 
time  instead  of  position  as  done  previously.  Here  it  can  be 
seen  that  the  tip  of  the  tooth  follows  the  static  deflection 
curve  for  each  of  the  speeds.  This  is  again  due  to  the  slow 
and  gradual  engagement  of  the  load. 

(4.4.4)  Dynamic  Deflections:  Rim  Included 

Including  the  rim  adds  flexibility  to  the  system  as 
already  mentioned.  The  tooth  tip  deflection  history  for  the 
impact  loading  case,  shown  in  Figure  4-15,  illustrates  this 
fact.  As  before,  the  amplitude  and  frequency  of  vibration  are 
the  same  for  each  of  the  four  speeds.  However,  compared  to  the 
beam  constraint  case  of  Figure  4-10,  the  amplitude  is  con- 
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siderably  larger  and  the  frequency  slower. 

When  the  finite  engagement  rise  time  loading  case  is  used 
the  same  general  results  are  noted  as  before.  That  is,  as  the 
speed  of  the  moving  load  increases,  the  rise  time  approaches 
impact  conditions  for  V*=0.01  (see  Figure  4-16). 

Application  of  the  Wallace-Seireg  load  history  equations 
to  the  rim  constraint  case  produces  results  which  behave 
exactly  as  before.  In  Figure  4-17  the  tip  of  the  tooth 
deflects  in  proportion  to  the  magnitude  of  the  applied  load. 
The  oscillations  about  the  static  curve  are  evident  but  do  not 
contribute  significantly  to  the  overall  response. 


5.  DETERMINING  THE  EFFECTS  OF  INERTIA  ON  THE 
DYNAMIC  RESPONSE  OF  MESHING  GEAR  TEETH 

In  most  theoretical  models  developed  to  determine  the 
dynamic  response  of  meshing  gear  teeth  [1-5],  the  mass 
(inertia)  of  the  tooth  is  neglected  in  the  analysis.  This 
simplification  is  based  on  the  claim  that  the  mass  of  the  tooth 
is  small  compared  to  the  mass  of  the  hub. 

To  determine  the  effect  of  the  inertia  of  the  tooth  on  the 
dynamic  response  of  a  meshing  gear  pair,  a  simplified  model  of 
two  cantilever  beams  attached  to  foundation  masses  is  used. 

The  system  analyzed  is  illustrated  in  Figure  5-1. 


Two  cantilever  beams  with  identical  geometric  and  material  pro¬ 
perties  are  rigidly  fixed  to  two  foundation  masses.  Ml  and  M2. 
In  this  analysis,  the  masses  of  the  foundations  are  defined  to 
be  the  same.  The  two  beams  are  held  together  by  opposing  for¬ 
ces,  P,  acting  on  the  masses  Ml  and  M2.  A  massless  ball  acts 
as  the  contact  point  between  the  beams  and  moves  along  the 


84 


v 


ft 


■  *1 

*1* 


I 


t 


beams  at  a  prescribed  speed.  As  the  contact  point  moves  from 
position  1  to  position  2  at  a  velocity  V,  the  change  in  stiff¬ 
ness  of  the  beams  causes  the  masses  to  oscillate  in  the  direc¬ 
tion  of  P  at  the  system  resonant  frequency.  To  simplify  the 
problem,  movement  only  in  the  direction  of  P  is  allowed. 

The  dynamic  response  of  the  meshing  cantilever  beams  is 
determined  for  two  loading  cases;  one  where  the  beams  are 
assumed  massless,  the  other  where  the  masses  of  the  beams  are 
included.  The  system  is  analyzed  using  constant  and  variable 
speed  moving  loads  of  constant  magnitude.  Both  impact  and 
smooth  load  engagement  responses  are  examined  by  changing  the 
initial  conditions  of  the  system.  These  loading  cases  are  ana¬ 
lyzed  using  two  values  for  the  foundation  masses  of  1.0  and 
1.0E-4  lbs. 


(5.1)  Analysis  Using  Massless  Beams 

Figure  5-2  shows  the  system  used  when  the  beams  are 
assumed  massless.  As  the  contact  point,  represented  by  a 


Figure  5-2:  System  parameters 
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massless  ball,  moves  along  the  beams,  the  deflection  of  the 
beams  will  change  due  to  the  variation  in  stiffness.  The 
stiffness  of  each  beam  varies  with  the  local  coordinate,  E,  . 
Since  the  beams  are  massless,  they  do  not  affect  the  system 
resonant  frequency,  but  follow  the  oscillations  of  Ml  and  M2 
exactly.  To  determine  the  dynamic  response  for  the  massless 
beam  configuration,  the  differential  equations  of  motion  are 
written  and  solved  for  this  system. 

(5.1.1)  Equations  of  Motion 

The  equations  of  motion  for  this  system  are  determined  by 
first  considering  each  beam-mass  configuration  as  a  free  body 
(see  Figure  5-3).  Writing  Newton's  Second  Law  as  the  sum  of 


Figure  5-3:  Free  body  diagrams 


the  forces  acting  on  each  body  we  have; 

Ml  XI  -  F  -  p  body  1  (5-1) 

M2  X2  =  -F  +  P  body  3  (5-2) 


(Here  it  is  assumed  that  XI  and  X2  define  the  positive  displa¬ 
cement  direction).  Dividing,  through  by  the  coefficient  of  the 


second  derivative  term,  and  subtracting  (5-1)  from  (5-2) 
yields; 

X2  -  XI  *  +  P(^+^)  C5-3) 

To  further  simplify  equation  (5-3),  the  right  hand  side  is  com¬ 
bined  and  a  common  denominator  is  determined.  This  gives; 

••  ••  M1+M2 

X2  -  XI  =  (P-F)£gjgp  (5-4) 

where; 

1  1  M1M2  ,1  1  .  M1+M2 
Ml  M2  M1M1  m  M2'  “  Ml  M 2 

Making  simple  substitutions,  equation  (5-4)  can  then  be  written 
as; 

MX  *  P-F  (5-5) 


where; 


and; 


M  = 


M1M2 

M1+M2 


X  *  X2-X1 


From  beam  theory,  the  static  deflections  of  each  beam  at 
the  point  of  contact  of  the  load,  are  given  by; 

F£l3  F523 

S1  -TB-  '  S2  =  W  (5-6> 

as  illustrated  in  Figure  5-4.  Equations  (5-6)  are  written  in 
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Figure  5-4:  Beam  deflection  configuration 


terms  of  the  beam  stiffness  as; 


(5-7) 


where; 

K1  =  S  ;  K2  =  ~ 

5l3  C23 

represent  the  stiffnesses  of  the  beams. 

As  the  load  moves  along  between  the  beams,  the  varying 
stiffness  causes  Ml  and  M2  to  oscillate  at  the  resonant  fre¬ 
quency  of  the  system.  To  insure  constant  contact  between  the 
beams  while  the  masses  are  vibrating,  the  following  rela¬ 
tionship,  termed  the  constraint  equation,  must  be  satisfied; 


X2  -  XI  «  61  +  62 


(5-8) 


where  61  and  62  assume  orientations  as  shown  in  Figure  5-4. 
Substituting  equations  (5-7)  into  (5-8)  gives  the  constraint 
equation  in  terms  of  the  beam  stiffnesses  as; 
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*2  -  X1  -  KT  +  &  (5-9> 

Combining  the  right  hand  terms  in  equation  (5-9), 

X2  -  XI  =  F(^~£p) 

multiplying  through  by 

v  -  K1K2 
K1+K2 

and  letting  X=X2-X1,  yields  a  familiar  form  of  the  constraint 
equation,  describing  the  deflection  of  a  spring; 

F  -  KX  (5-10) 

Inserting  F  from  equation  (5-10)  into  equation  (5-5)  gives  the 
equation  of  motion  for  the  massless  beam  system; 

MX  +  KX  -  P  (5-11) 


where; 


M  ••  •• 

X  ■  X2  -  XI 
X  *  X2  -  XI 

„  _  M1M2 
”  *  M1+M2 

_  ,  K1K2 
K  K1+K2 

Initially,  the  system  deflection  is  determined  by  assuming  that 
static  equilibrium  is  satisfied.  Thus,  at  time  equal  to  zero, 
the  initial  deflection  is  the  applied  load  divided  by  the  total 
stiffness  (the  combined  deflections  of  the  beams); 
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X(0)  =  X2 ( 0 ) -XI (0 )  =  P/K 


With  the  initial  deflections,  equation  (5-11)  is  solved 
using  a  fourth  order  Runge-Kutta  integration  algorithm.  Since 
the  Runge-Kutta  algorithm  used  is  designed  for  systems  of  first 
order  equations,  the  second  order  differential  equation  (5-11) 
must  be  converted  to  first  order.  This  is  done  by  defining; 

Y1  *  X  and  Y2  *=  X 

Substituting  these  relations  into  (5-11),  we  then  have; 

MYl  +  KY2  =  P  (5-12) 


where  at  t=0.0. 


Y1(0)  =  0.0 
Y2  (0 )  =  PA 

Given  these  initial  values,  the  relations; 

Y1  *  (P-KY2 )/M 
Y2  -  Y1  »  X 

are  used  as  input  for  Runge-Kutta  and  integrated  to  determine 
the  displacement,  X-X2-X1  during  the  load  cycle.  (See  Appendix 
7  for  solution  algorithm). 

(5.1.2)  Static  Analysis 

Before  any  dynamic  analysis  is  performed,  a  static 
deflection  test  is  done  to  provide  a  reference  for  dynamic 
deflection  comparisons.  Solving  the  relationship; 
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X2  -  XI  -  P/K 


at  different  load  positions  produced  the  plots  shown  in  Figure 
5-5.  The  three  blocks  represent  the  components  comprising  the 
constraint  equation  (5-8).  Deflections  are  determined  between 
load  positions  1  and  2  illustrated  in  Figure  5-6.  The  load 
position  is  measured  relative  to  beam  1  as  labelled  on  the 
abscissa  of  Figure  5-5. 

(5.1.3)  Dynamic  Response  Results 

The  dynamic  response  of  the  massless  beam  system  is  deter¬ 
mined  for  two  loading  cases;  a  constant  speed  1000  lb  load,  and 
a  variable  speed  1000  lb  load.  For  each  of  these  loading 
cases,  two  sets  of  initial  conditions  are  considered;  those 
defined  for  static  equilibrium  in  equation  (5-11),  and  another 
set  where  the  initial  deflection,  X2  -  XI,  is  equal  to  zero. 
The  second  of  these  initial  conditions  will  cause  the  load  P  to 
be  experienced  as  an  impact  load  since  the  beams  initially  will 
have  no  deflection  and  will  attempt  to  return  to  static 
equilibrium. 

For  loads  moving  with  constant  speed,  speeds  of;  1.0,  5.0, 
10.0,  20.0,  and  40.0  inches  per  second  are  used.  This  range  of 
speeds  is  chosen  in  order  to  bracket  the  cycle  period  asso¬ 
ciated  with  the  fundamental  frequency  of  the  system.  The  fre¬ 
quency  of  vibration  of  the  system  is  constantly  changing  with 
the  position  of  the  load.  However,  a  representative  fundamen¬ 
tal  frequency  is  calculated  when  the  load  is  halfway  through 
the  load  cycle.  At  this  position  the  beam  stiffnesses  are 
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equal.  Then,  we  simply  have  a  spring  mass  system  whose 
undamped  natural  frequency  is  defined  as; 


0) 


(cycle/sec) 


With  beam  dimensions  the  same  as  those  given  in  Figure  5-6  and 
the  masses  Ml  and  M2  equal  to  1.0,  the  period  of  oscillation  is 
approximately  0.03  seconds.  Thus  for  a  constant  velocity  of 
7.0  in/sec,  approximately  two  oscillations  will  occur  during  the 
load  cycle. 

When  a  variable  speed  moving  load  is  used,  the  speed  ini¬ 
tially  is  zero  and  increases  linearly  to  1.0,  5.0,  10.0,  20.0, 

and  40.0  in/sec  at  the  end  of  the  cycle, 
i 

(5. 1.3.1)  Constant  Speed  Moving  Loads 

For  slowly  moving  loads,  relative  to  the  fundamental 

I 

i 

>  period  of  the  system,  one  would  expect  near  static  deflections. 

i 

^  This  is  the  case  for  a  speed  of  1.0  in/sec  as  shown  in  Figure 

5-7.  The  figure  gives  the  components  of  the  constraint 
equation  resulting  from  an  initial  deflection  equal  to  the 
ratio  of  the  load  and  stiffness.  The  only  deviations  from  the 
static  deflection  curve  are  caused  by  small  oscillations  at 
the  resonant  frequency.  As  the  speed  of  the  moving  load 
increases,  one  would  expect  an  j ncrease  in  the  amplitude  of 
oscillation  as  the  period  of  the  load  cycle  approaches  the 
resonant  period.  Increasing  to  5.0  in/sec  produced  larger 
deviations  from  static  due  to  increased  oscillation  amplitudes 
(Figure  5-8).  This  trend  continues  until  the  deflection  is 
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Figure  5-7 
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maximum  for  a  speed  of  20.0  in/sec.  (See  Figure  5-9  and  5-10). 
In  both  Figures  5-9  and  5-10 ,  the  value  of  X2-X1  goes  negative. 
This  means  that  the  beams  have  separated.  Notice  also  that  the 
system  at  20  in/sec,  does  not  even  begin  to  approach  the  ini¬ 
tial  conditions  at  the  end  of  the  cycle,  as  for  the  slower 
speeds.  A  further  increase  in  speed  to  40.0  in/sec,  shows  that 
although  the  beams  themselves  deflect  appreciably  ( DELTA1  and 
OELTA2 ) ,  the  masses  themselves  are  displaced  very  little 
because  the  load  cycle  is  much  shorter  than  the  resonant  period 
(see  Figure  5-11). 

In  each  case,  the  system  strives  toward  the  static  deflec¬ 
tion  position.  But  due  to  the  inertia  of  the  foundation 
masses,  this  position  may  or  may  not  be  maintained  depending  on 
the  speed  of  the  moving  load. 

As  discovered  in  the  finite  element  analysis  (see  section 
4.4.3),  the  type  of  load  engagement  significantly  affects  the 
response  of  the  undamped  gear  tooth.  This  is  also  the  case  for 
the  meshing  beams  configuration.  By  changing  the  initial  con¬ 
ditions  of  X-X2-X1  to  zero,  the  system  reacts  to  restore  itself 
to  static  equilibrium.  Since  the  system  is  undamped,  a  high 
amplitude  oscillation  is  set  up  due  to  the  rapid  movement  of 
the  foundation  masses.  This  is  best  illustrated  by  examining 
Figure  5-12  where  the  components  of  the  constraint  equation 
are  plotted  for  a  velocity  of  1.0  in/sec.  As  before,  the 
system  oscillates  about  the  static  deflection  position,  but 
with  very  large  amplitude.  Note  also,  that  separation  occurs 
during  each  oscillation.  The  same  phenomenon  occurs  for 
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speeds  of  5.0  and  10  in/sec  (Figure  5-13  and  5-14).  However, 
with  a  further  increase  in  speed,  the  system  does  not  have  suf¬ 
ficient  time  to  react  and  the  oscillations  become  less  signifi¬ 
cant  (see  Figure  5-15  and  5-16).  One  should  not  be  led  to 
believe  that  the  separation  occuring  for  these  initial  con¬ 
ditions  is  caused  by  the  characteristics  of  the  system.  It  is 
caused,  however,  by  the  application  of  a  1000  lb  load  to  a  beam 
which  in  the  physical  sense  could  not  support  it.  However, 
since  the  system  behaves  linearly,  the  characteristics  of  the 
response  are  the  same  for  whatever  load  is  applied,  and  the 
1000  lb  load  is  used  simply  to  exaggerate  that  response. 

(5. 1.3. 2)  Variable  Speed  Moving  Loads 

Introducing  a  variable  speed  moving  load,  versus  constant 
speed,  has  little  effect  on  the  general  behavior  of  the  meshing 
massless  cantilever  beam  system.  The  same  conclusions  can  be 
drawn  concerning  the  dynamic  resonse  trends  due  to  increased 
moving  load  speed.  As  a  reference,  the  dynamic  response  curves 
for  both  sets  of  initial  conditions  and  the  five  different 
speeds,  plotted  as  a  function  of  position,  are  included  in 
Appendix  4. 

(5.2)  Analysis  Including  Beam  Inertia 

Figure  5-17  illustrates  the  physical  system  used  to  deter¬ 
mine  the  effects  of  including  the  inertia  of  the  beams  on  the 
dynamic  response  due  to  moving  loads.  This  system  is  exactly 
the  same  as  the  one  used  in  Section  5.1  except  the  mass  of  the 
beams  is  included.  First,  the  equations  of  motion  for  this 
system  are  developed. 
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Figure  5-17:  System  parameters 


(5.2.1)  Equations  of  Motion 

The  equations  of  motion  for  this  system  are  somewhat  more 
difficult  than  those  developed  for  the  massless  beam  con¬ 
figuration.  Instead  of  writing  down  the  differential  equation 
directly,  a  form  of  Lagange's  Equation  is  used  which  takes  into 
account  the  added  mass  of  the  beams. 

The  Lagrange  Equation  of  motion,  utilizing  Lagrange 
multipliers  is; 

d  f  9L  v  9L  _  n  ?  x 
dt  9qi  9qi  Qi  k=1/kaki 

for:  i  -  1,2,. ..,n  (5-13) 


where; 


L  -  T  -  0  »  Lagrangian 
T  *  System  kinetic  energy 
U  *  System  potential  energy 
qi  -  Generalized  coordinates 
Qi  »  Sum  of  non-potential  forces 
A]t  »  Lagrange  Multipliers 


aki 


9  f (gl ,  g2 , « » » , gn ) 
3qi 


first  partial 


derivative  of  the  constraint  equation 
with  respect  to  the  generalized  coordinates 


The  last  term  on  the  right  hand  side  of  equation  (5-13)  repre¬ 
sents  the  constraint  forces  which  makes  it  possible  to  regard 
all  generalized  coordinates,  qi,  as  independent. 

Four  generalized  coordinates  are  used  to  describe  the 
system  position,  velocity,  and  acceleration.  They  are; 

XI  -  describing  Ml 
X2  -  describing  M2 
ql  -  describing  beam  2  relative  to 
foundation  mass  Ml 
q2  -  describing  beam  2  relative  to 
foundation  mass  M2 


From  the  definition  of  generalized  coordinates,  these  four 
quantities  are  assumed  independent  of  each  other,  and  when  used 
together  they  describe  the  state  of  the  dynamical  system  at  any 


The  static  deflection  of  a  loaded  cantilever  beam  is  given 
by  the  applied  load  divided  by  the  stiffness.  When  considering 
the  dynamic  deflection  of  a  beam  which  has  mass,  this  rela¬ 
tionship  is  not  valid.  Instead,  the  superposition  of  the 
natural  vibrational  modes  is  used  to  describe  the  shape  of  the 
deformation,  and  when  multipled  by  the  generalized  coordinate, 
qi,  the  actual  deflection  is  obtained. 

The  mode  shapes  are  an  infinite  set  of  eigenvectors 
derived  from  the  differential  equation  of  the  cantilever  beam. 
To  exactly  duplicate  the  deflection  obtained  from  elasticity 
theory,  an  infinite  number  of  modes  must  be  used.  In  practice, 
however,  only  the  first  few  contribute  significantly  to  the 
overall  deflection,  such  as  those  illustrated  in  Figure  5-18. 


Figure  5-18:  Natural  Modes 


In  this  analysis/  only  the  first  node  shape  is  used/  the 
natural  mode  described  by  the  equation; 

<P(C)  *  Cl  (sinBL-sinhBL)  ( sin 6£-sinh3£ ) 

+  ( cos  8L+  cosh BL)  (cosB£-coshB£)  ]  (5-14) 

where  the  parameters  are  defined  in  Figure  5-18. 

We  are  now  ready  to  develop  the  equations  of  motion  for 
the  mass-beam  system  of  Figure  5-17.  The  kinetic  energy  of  the 
system  is  determined  by  considering  the  velocity  of  the  foun¬ 
dations  Ml  and  M 2,  and  the  vibration  of  the  beams  with  respect 
to  the  foundation  masses.  By  defining  the  beam  deflections  to 
be  positive  as  shown  in  Figure  5-19/  the  kinetic  energy  of  the 

X  2 - - 


P 


- P 

Figure  5-19:  Beam  deflection 
configuration 


system  can  be  written  as; 

T  =  JmI  XI2  +  4«2  X22  + 


1  /■kl.  . 

2**-  '  2*»“  ““  •  j  P  I  41!  qj_)  d^i 

J0 


1  l*2* 

JP  r  (X2  -  4>2  <l2)2d5] 
Jn 


(5-15) 


no 


The  potential  energy  consists  only  of  the  strain  energy  due  to 
the  beam  deflections; 


u 


=  |ei 


r 


(4>i"  qi)^!  + 


{ 


L2 


(4>2"  <32>  d^2 


(5-16) 


Where  4)”  is  the  second  derivative  of  the  mode  shape  with 
respect  to  the  length  variable,  £  .  In  equations  (5-15)  and 
(5-16)  the  mass  per  unit  length,  p  ,  the  elastic  modulus,  B, 
and  the  moment  of  inertia,  I,  are  constant  over  the  length  of 
the  beams  and  can  be  left  outside  the  integrals.  Substracting 
(5-16)  from  (5-15)  to  form  the  Lagrangian  we  have; 


L  «  T  -  0  «  ^ Ml  XI2  +  ^M2  X22 


1  #  O  *1  a 

JP  j  (Xl  +  <}>!  Sirdq  +  J  Pj  (X2  -  $2  q2) 

Jo  Jo 

-  Iei  -  5*1 


(5-17) 


For  the  meshing  cantilever  beam  system  there  are  four 
unknown  parameters,  XI,  X2,  ql,  q2  and  their  derivatives,  each 
dependent  on  time,  describing  the  dynamical  system.  At  any 
time,  t,  the  position  of  the  system  can  be  described  by  the 
constraint  equation; 

X2-X1  -  -  4>2q2  =  f  (Xl,X2,ql,q2)  -  0  (5-18) 

If  the  constraint  equation  is  not  satisfied,  separation  occurs 
between  the  beams  at  the  point  of  contact. 

A  Lagrange  Equation  of  motion  is  written  for  each  degree 


of  freedom  of  the  system  corresponding  to  the  four  generalized 


coordinates.  This  gives  four  equations  relating  the  genera- 
lized  coordinates  and  X  .  Including  the  constraint  equation  or 
its  derivative  gives  a  fifth  equation  necessary  to  determine 

In  terms  of  XI,  the  equation  of  motion,  term  by  term, 
becomes ; 

=  M1X1  +  p  f  (XI  + 


4r(4£-)  =  Ml  XI 
dt  3X1 


+  P 


n 


(XI  + 


*iqi 


)dC- 


(5-19a) 


and; 


-  0.0  (5-19b) 

On  the  right  hand  side,  differentiation  of  the  constraint 
equation  with  respect  to  XI  yields; 

$r--1.0  (5-190 

Combining  equations  (5-19a)  through  (5-19c)  along  with  the  non¬ 
constraint  force,  P,  acting  on  Ml,  in  the  form  of  equation 
(5-13)  we  have; 


(Ml  +  p 


r 


The  term  p 


/'Ll 

Jo 


defined  as  MB1 . 
coordinate; 


1 

dCx)Xl  +  p 


*1>«1 


P-X 


(5-20) 


is  simply  the  total  mass  of  beam  1  which  is 
Similarly,  with  X2  as  the  generalized 
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(5-21) 


(M2  +  MB2)X2 


P  +  X 


Performing  the  required  differentiation  of  the  Laqrangian  and 
constraint  equation  with  respect  to  ql  yields; 


9L 

3qi 


fL1  *  2 

pj^  XI  + 


qi)d5] 


*  PJ?<*1  il  +  *l2  i'l,d51 

3f  -  * 

^T‘  *1 


(5-22a) 

(5-22b) 

(5-22c) 


Combining  equations  (5-22a)  through  (5-22c),  and  noting  that  no 
external  (non-constraint)  forces  act  on  the  beams,  the  equation 
of  motion  derived  with  respect  to  ql  becomes; 


/■LI  /Ll  j  Ll 

p  Jo  (^XDd^+pJ  (<(>!%)  d^1+EIJ0  l 2  <Jl)d^i  ■  "♦ix 


rLl 

(♦.XDdC^+P  I 
io  1  1  Jo 

Likewise  for  q2; 


(5-23) 


-P  JLU2X2)d?2+P  JLU22q2)d£2+EI  fL*<P$2  q2)dC2  =  - 


<j>2X 


(5-24) 


The  fifth  equation  of  motion  is  determined  by  differentiating 
the  constraint  equation  twice  with  respect  to  time.  This  gives 
an  equation  relating  the  second  derivatives  of  the  generalized 
coordinates  with  respect  to  time.  Differentiating  equation 
(5-18)  once  yields; 


1)3 


X2-Xl-$i  qi”^2  92- *2  92*  0 


After  a  second  differentiation  and  rearrangement; 

)C2-Xl-$i  9l“  $2  92  *  ^1  9l+^2  92+^^l  91+^2  92  ^  (5-25) 

Equations  (5-20),  (5-21),  (5-23),  (5-24)  and  (5-25)  are  the 

equations  of  motion  describing  the  dynamical  system  of  Figure 
5-16. 

These  equations  can  be  greatly  simplified  by  making  the 
following  substitutions: 

PHI1  »  ;  PHI  2  =  <p2 

rLl  rL2 

IPHI1  =  p  ?  IPHI2  =  p  <P2d^2 

-'o  Jo 

{•LI  *L2 

p  (j»12dC1  *  1  ;  P  $2  d?2  =  1 

J0  J0 

rhl  fL2 

I2PHI1  -  El  J  dC!  r-  I2PHI2  *  El  *22&Z2 

F(q1»q2»^i»,i,2)  “  $1  9X  ♦  $2  92  +  2^1  9i  +  ^2  ^2) 

(see  Appendix  5  for  developments  concerning  these 
simplifications).  The  equations  can  then  be  written  in  matrix 
from  in  terms  of  the  second  derivatives  and  as; 


lAKX)  -  { B> 


(5-26) 


M  i  « 


■"& VjjV-t  .4  i*K  ;V  *  V 


where; 


[A]  * 


( M1+MB1 ) 
0 

IPHI1 

0 

-1 


0 

IPHI1 

0 

1 

M2+MB2 

0 

-IPHI2 

-1 

0 

1 

0 

PH  11 

-IPHI2 

0 

1 

PBI2 

1 

-PHI1 

-PHI2 

0 

(X)  = 


and  {B}  * 


-P 
P 

-I2PHI*ql 
t -I2PHI2*q2 
F(ql,q2,4>l,<t>2)l 


The  only  differences  between  equations  (5-26)  and  the  equation 
of  motion  for  the  massless  beam  system  (equation  (5-11))  are 
the  inertial  terms.  These  are;  X1*IPHI1,  X2*IPHI2,  and  the 
terms  involving  ql  and  q2.  By  eliminating  the  inertial  terms 
from  matrix  [A],  the  resulting  equations  describe  the  massless 
beam  system.  Appendix  6  includes  the  analysis  and  results  from 
this  test. 

The  initial  conditions  governing  the  system  described  by 
equations  (5-26)  are  chosen  such  that  the  massless  beam  system 
is  emulated.  When  chosing  the  initial  conditions ,  there  are 
two  sets  of  parameters  which  must  be  considered.  The  first  set 
consists  of  the  four  position  variables  XI,  X2 ,  ql  and  q2.  In 
section  (5.1.1)  it  is  stated  that  the  beams  assume  a  static 
orientation  at  the  beginning  of  the  load  cycle.  Thus,  by 
chosing  ql  and  q2  as; 


ql  (0) 


P*PHI1 
*  I2PHI1 


q2  (0) 


P*PHI2 

I2PHli 


(5-27) 
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and  XI (0)  equal  to  zero  as  before,  the  deflection  X2(0)  is 
determined  from  the  constraint  equation  as; 

X2 (0 )  -  ql (0) PHI1  +  ql (0) PHI2  (5-28) 

This  same  reasoning  can  be  applied  to  the  first  derivatives  of 
the  generalized  coordinates,  XI,  X2 ,  ql  and  q2.  Previously, 
the  initial  velocities  of  the  foundation  masses  were  determined 
such  that  X2-X1*0  at  time  equal  zero.  For  this  case,  the  same 
effect  is  obtained  by  forcing; 

X1(0)  *  0.0 

X2(0)  *  0.0  (5-29) 

ql (0)  -  0.0 

Then,  using  the  first  derivative  of  the  constraint  equation  the 
value  of  q2(0)  is  determined  from; 

X2  -  XI  -  q j*PHIl  +  q2*DlPHIl  +  q2*PHI2  +  q2*DlPHI2 

where  at  time  equal  to  zero,  all  terms  are  zero  except  for 
those  in  equation  (5-30); 

q2(0)  -  ( -qj* DlPHIl  -  q2* D1PHI2)/PHI2  (5-30) 

where  DlPHIl,  D1PHI2  and  PHI2  are  evaluated  at  the  beginning  of 
the  load  cycle  and  ql  and  q2  are  taken  from  equations  (5-27). 

With  the  initial  conditions  described  by  equations  (5-27) 
through  (5-30),  the  dynamic  response  of  the  system  described  in 
equations  (5-26)  is  determined. 

For  each  time  step  during  the  load  cycle,  the  system  of 
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equations  (5-26)  is  inverted  to  determine  the  column  matrix 
{ X>  *  { XI  X2  ql  q2  X}  T  in  terms  of  the  current  values  of  { B} 
on  the  right  hand  side.  At  the  beginning  of  the  load  cycle  {X> 
is  determined  for  those  initial  conditions  set  forth  pre¬ 
viously,  along  with  the  right  hand  side  (B>.  The  values  of  {  X} 
are  then  used  with  the  first  derivatives  of  the  generalized 
coordinates  (initially  equations  (5-29)  and  (5-30))  as  input 
for  a  Runge-Kutta  integration  routine  which  integrates  for  the 
desired  solution  parameters  XI,  X2,  ql  and  q2.  It  also  deter¬ 
mines  the  new  XI,  X2,  ql  and  q2  which  are  in  turn  used  for  the 
next  integration  step.  The  algorithm  containing  the 

Runge-Kutta  ( RKGS )  and  matrix  inversion  (SIMQ)  subroutines 
along  with  the  parameters  describing  equations  (5-26)  is 
included  in  Appendix  7.  Therein  lies  a  detailed  description  of 
the  programming  steps  for  the  solution  of  equations  (5-26) 
using  the  aforementioned  initial  conditions. 

(5.2.2)  Dynamic  Response  Results:  Foundation  Mass  *  1.0  lbs 
As  stated  in  section  (5.2.1),  only  the  first  material  mode 
of  vibration  is  included  when  determining  the  dynamic  response 
of  the  meshing  cantilever  beams  where  the  inertia  of  the  beams 
is  included.  One  would  then  think  that  the  system  whose  dyna¬ 
mic  response  is  composed  of  a  single  mode  shape  would  be 
somewhat  stiffer  than  the  same  configuration  where  the  beam 
deflection  is  determined  from  elasticity  theory.  However,  com¬ 
paring  Figures  5-7  to  5-11,  with  A6-1  to  A6-5,  it  is  seen  that 
the  difference  in  deflections  is  negligible. 
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In  this  analysis,  with  the  beam  inertia  included,  the 
dynamic  response  is  compared  to  the  massless  beam  problem.  At 
slow  moving,  constant  speed  loads,  the  dynamic  response  con¬ 
sists  of  quasi-static  deflection  with  the  moving  load  response 
superimposed  over  it  (see  Figure  5-20).  (Results  for  moving 
load  speed  of  1.0  in/sec  were  not  obtainable  due  to  lack  of 
convergence  in  the  Runge-Kutta  integration  routine). 
Continuing  to  increase  the  moving  load  speed  causes  an  increase 
in  the  response  as  the  resonant  frequency  is  further  excited  as 
shown  in  Figure  5-21  and  5-22.  Increasing  the  moving  load 
speed  above  the  system  resonant  frequency  produces  a  much 
smoother  response  as  illustrated  in  Figure  5-23.  Now,  examina¬ 
tion  of  Figures  5-20  through  5-23  compared  to  Figures  5-7 
through  5-11  shows  that  for  values  of  the  foundation  masses  of 
1.0  lbs  and  moving  load  speeds  of  5.0,  10.0,  20.0  and  40.0 
in/sec,  the  responses  of  the  two  cantilever  beam  systems  are 
essentially  the  same.  The  same  conclusions  are  drawn  when  exa¬ 
mining  Figures  5-24  through  5-27  as  compared  to  Figures  5-13 
through  5-16.  Here,  the  initial  conditions  are  changed  such 
that  the  beams  experience  an  impact  load  at  the  beginning  of 
the  load  cycle.  From  the  figures  it  can  be  seen  that  the 
responses  of  the  two  systems  are  again,  very  much  the  same. 

As  done  for  the  massless  beam  configuration,  the  speed  of 
the  moving  loads  is  allowed  to  vary  linearly  from  zero  to  a 
maximum  value  during  the  load  cycle.  This  analysis  serves  as  a 
useful  check  since  the  initial  conditions  are  easily  defined 
due  to  a  stationary  load.  The  results  from  this  analysis  are 
contained  in  Appendix  4. 
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(5.2.3)  Dynamic  Response  Results: 

Foundation  Mass  «  1.0E-4  lbs. 

In  order  to  substantiate  the  conclusion  that  the  dynamic 
response  is  not  dependent  on  the  mass  of  the  beam,  both  the 
massless  beam  and  inertial  beam  configurations  are  analyzed 
again  with  a  decrease  in  the  difference  between  the  values  of 
the  beam  and  foundation  masses.  Of  course  the  mass  of  the  beam 
is  not  included  in  the  massless  beam  analysis.  Previously  the 
mass  of  the  foundation  was  1.0  lbs  and  the  total  mass  of  the 
beam  was  approximately  3.6E-6  lbs.  In  this  analysis,  the 
foundation  mass  is  decreased  to  1.0E-4  lbs  while  the  mass  of 
the  beam  remains  the  same.  Since  the  overall  mass  of  the 
system  is  greatly  decreased,  the  corresponding  resonant  fre¬ 
quency  is  much  larger.  Using  an  approximation  of  the  resonant 
frequency  using  a  single  mass  and  the  load  at  the  midpoint  of 
the  beam  yields  an  approximate  resonant  frequency  of; 

w  *71  ^  -  3500  cycle/sec. 

and  a  cycle  period  of  approximately  2.875E-4  sec.  Therefore, 
in  order  to  provide  an  adequate  range  of  moving  load  speeds 
such  that  the  system  resonance  is  bracketed,  speeds  of  100., 
500.,  1000.,  1500.  and  2000.  in/sec  are  used.  (Approximately 
one  oscillation  cycle  occurs  at  1700  in/sec).  To  facilitate 
comparisons  of  respective  systems,  the  results  are  presented  in 
a  parallel  fashion.  Figures  5-28  and  5-29  show  the  respective 
responses  for  a  moving  load  speed  of  500  in/sec.  As  seen  from 
the  plots,  the  responses  are  virtually  the  same.  The  results 


are  similar  for  the  other  moving  load  speeds  (see  Figure  5-30 
through  5-35). 
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6 .  CONCLUSIONS 


The  objectives  of  this  analysis  were  twofold.  First,  the 
dynamic  response  of  a  single  spur  gear  tooth  subject  to  moving 
loads  was  studied  to  determine  the  effect  of  the  speed  of  move¬ 
ment  of  the  load.  A  spur  gear  tooth,  modelled  using  finite 
element  techniques,  was  used  in  the  analysis.  From  this  analy¬ 
sis  it  was  found  that  the  dynamic  response  of  a  single  gear 
tooth  is  not  dependent  on  the  speed  of  the  moving  load,  but 
rather  on  the  type  of  load  engagement  experienced  at  the 
beginning  of  the  load  cycle.  Including  the  rim  in  the  analysis 
added  flexibility  to  the  system  but  did  not  change  the  overall 
response. 

The  second  objective  was  to  determine  whether  or  not  the 
mass  (inertial  forces)  of  the  tooth  can  be  neglected  when  small 
compared  to  the  mass  of  the  gear  hub.  A  simplified  analysis 
using  meshing  cantilever  beams  was  used.  For  the  range  of 
speeds  tested,  it  was  found  that  the  inertia  forces  of  the 
tooth  are  small  and  therefore,  the  mass  of  the  tooth  can  be 
neglected  when  determining  the  dynamic  response  of  a  meshing 
gear  system. 
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APPENDIX  1 

1-  Finite  Element  Model 
Generation  Program 

2.  Moving  Load  Generation  Program 
Using  Wallace-Seireg  Load  History 
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SPUR  GEARTOOTH  PROFILE  GENERATOR 


* 

* 

* 


C  * 

C  * 

c  * 


:  v: 


t'* 


Uf 


\  i 


IT 


c  * 

C  *  WRITTEN  BY  LYLE  W.  SHUEY 

C  * 

C  *  THIS  PROGRAM  GENERATES  A  FINITE  ELEMENT 

C  *  MESH  FOR  EITHER  A  HIGH  CONTACT  RATIO 
C  *  GEAR  OR  A  LOW  CONTACT  RATIO  GEAR,  IN- 

C  *  CLUDING  THE  FOUNDATION  AND  RIM  IN  THE 

C  *  MODEL. 

C  *  A  MINIMUN  OF  INPUT  PARAMETERS  ARE  NEEDED 


c 

* 

TO  DESCRIBE  THE  GEAR  GEOMETRY. 

c 

* 

THEY  ARE; 

c 

* 

PRESSURE  ANGLE 

:PANG 

c 

* 

PITCH  RADIUS 

:RP 

c 

* 

ADDENDUM 

:  AD 

c 

* 

DEDENDUM 

:  DED 

c 

* 

CIRCULAR  PITCH 

:CIRP 

c 

* 

BACKLASH 

: BACKL 

c 

* 

FILLET  RADIUS 

:RF 

c 

* 

RIM  THICKNESS 

:RTH 

c 

* 

C  ****************************************** 

C 

c 

C******  the  program  reads  the  input  variables 

c******  USING  THE  PROCEEDING  STATEMENTS. 

C 

DIMENSION  DUMY (20) ,YPRF(20) ,XPRF(20) 
DIMENSION  XX(40) ,YY(40) ,XN(400) ,YN(400) 
DIMENSION  XANG(9,9) , XADA(9) ,XINC(9) 
PI-3.141592654 
READ (5,*)  PANG 
WRITE (6,*)  'PRESSURE  ANGLE: 

PANG— PANG*  PI/ 180 . 

READ(5, *)  RP 

WRITE (6,*)  ‘PITCH  RADIUS: 

READ (  5 ,  * )  AD 
WRITE (6,*)  'ADDENDUM: 

READ(5, *)  DED 
WRITE ( 6 , * )  • DEDENDUM : 

READ(5, *)  CIRP 
WRITE (6,*)  'CIRCULAR  PITCH: 

READ(5, *)  RF 

WRITE (6,*)  'FILLET  RADIUS: 

READ(5, *)  BACKL 
WRITE (6,*)  'BACKLASH: 

READ(5, *)  RTH 
WRITE (6,*)  'RIM  THICKNESS: 

RB— RP*COS (PANG) 

THP- ( (RP*RP-RB*RB)/ (RB*RB) ) ** . 5 
CIRTH— CIRP/2 . -BACKL 

WRITE (6,*)  'CIRTH:  CIRTH 

AL-THP-ATAN (THP) +CIRTH/ (RP*2 . ) 

DP— PI/CIRP 

WRITE (6,*)  'DIAMETRAL  PITCH:  '  ,DP 

RR-RP-DED 
RO-RP+AD 
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,  PANG 

,RP 
,  AD 
,  DED 
,CIRP 
,RF 

, BACKL 
,RTH 


Vi  ‘ 


OOOOI  ooooooooo 


*  THIS  SECTION  CALCULATES  POINTS  ON  THE  * 

*  INVOLUTE  SURFACE ,  FOUNDATION,  AND  RIM  * 

*  WHICH  WILL  LATER  BE  USED  TO  GENERATE  * 

*  NODAL  COORDINATES.  (BOTH  HCRG  AND  LCRG) * 


******  GENERATING  POINTS  ON  THE  INVOLUTE  PROFILE 


DO  5  1-1,11 
IA-1+ (1-1) 

HH— (RO-RB) /10. 

H— RO-HH* (1-1) 

ANG- ( (H*H-RB*RB) / (RB*RB) )  ** .  5 
X— RB*  (SIN  (ANG)  -ANG*COS  (ANG) ) 

Y-RB* (COS (ANG) +ANG*SIN (ANG) -1) 

ANGF-ATAN (SIN (ANG) /COS (ANG) ) -AL 

ANGF— ANGF*180 ./PI 

XP— RB*SIN (AL) +X*COS ( AL) +Y*SIN ( AL) 

YP— RB*COS (AL) -X*SIN (AL) +Y*COS ( AL) 

XX ( IA) — XP 
YY(IA)— YP 

WRITE ( 6 , * )  'XP: • ,XP, '  YP: • , YP, 'THETA: • , ANGF 
CONTINUE 

*  CHECK  IF  GEAR  IS  LCRG  OR  HCRG 


CHECK— RR*RR+2 . *RR*RF 
ITAN-2 

IF ( CHECK . LT . RB*RB)  ITAN-1 
IF(ITAN.NE.l)  GO  TO  15 
C 

C******  STRAIGHT  LINE  TANGENT  TO  INVOLUTE  (HCRG) 
C 


DD- ( -2 . *RR+ ( 4 . *RR*RR+4 . *2 . *RF*RR) ** . 5) /2 . 
D=RB- (RR+DD) 

DDD=RB-RR 
XTAN— ( D/DDD) 

XFILL- (DD/DDD) 

LTAN-INT (XTAN*8 . ) 

LFILL-INT (XFILL*8 . ) +1 
WRITE ( 6 , * )  'DD: 


•DJ 

DDD: 

'XTAN: 

•XFILL: 

•LFILL: 

•LTAN: 


' ,  DD 
'  •  D 
'  ,  DDL 
' , XTAN 
', XFILL 
' , LFILL 
' , LTAN 


WRITE ( 6 , * ) 

WRITE ( 6 , * ) 

WRITE ( 6 , * ) 

WRITE ( 6 , * ) 

WRITE (6,*) 

WRITE (6,*) 

DRI-D/LTAN 
DR— 0 . 0 

DO  10  I -1, LTAN 
IA-12+ (1-1) 

XP-XP-( ( (RB-DR) *SIN (AL) ) - ( (RB-DRI) *SIN(AL) ) ) 
YP-YP- ( ( (RB-DR) *COS (AL) ) - ( (RB-DRI) *COS (AL) ) ) 


DR-DRI 

DRI-DRI+D/LTAN 
XX(IA)« XP 
YY(IA) -YP 
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CONTINUE 


10 
C 

C*********  FILLET  RADIUS  SECTION  (HCRG) 

C 

BETA-ATAN ( (RR+RFJ/RF) 

BETA*  BETA/  LF I LL 
DO  11  1*1 , LFILL 

IA* ( 12+LTAN) + (1-1) 

X-RF-RF*COS (BETA* I) 

Y— RF*SIN  (  BETA*  I) 

XP- (RB-D) *SIN (AL) +X*COS ( AL) +Y*SIN ( AL) 

YP- (RB-D) *COS (AL) -X*SIN (AL) +Y*COS (AL) 
XX(IA)-XP 
YY(IA)»YP 
11  CONTINUE 

GO  TO  21 
C 

C******  GENERATING  POINTS  FOR  FILLET  RADIUS  (LCRG) 
C 

15  DD=RB-RR 

RF- ( ( DD*DD) + ( 2 . *DD*RR) ) / ( 2 . *RR) 

WRITE  (6,  *)  'NEW  FILLET  RADIUS  ',RF 
BETA=ATAN ( (RR+DDJ/RF) 

WRITE (  6 ,  * ) 'BETA: ' ,BETA*180./PI 

BETA* BETA/ 8 . 0 
DO  20  1*1,8 
IA*12+ (1-1) 

X*RF-RF  *  COS (BETA* I) 

Y*-RF*S IN ( BETA* I ) 

XP=RB*SIN (AL) +X*COS (AL) +Y*SIN (AL) 

YP=RB*COS (AL) -X*SIN (AL) +Y*COS (AL) 

XX ( IA) =XP 
YY(IA)*YP 

WRITE ( 6 , * )  'XP: • ,XP, 'YP: • ,YP 

20  CONTINUE 
C 

C******  GENERATING  POINTS  ON  OUTER  RIM  SURFACE 
C 

21  ADA*ASIN( (XP+CIRTHJ/RR) -ASIN(XP/RR) 
DADA-ADA/6. 

ADA1-ASIN (XP/RR) 

WRITE ( 6 , * ) 'ADA: ' ,ADA*180./PI 
DO  30  1*1,4 
IA*20+ (1-1) 

DUMY (I) *10 . 

XP*XP+RR* (SIN (ADA1+DADA) -SIN(ADAl) ) 
YP*YP-RR* (COS (ADA1) -COS (ADA1+DADA) ) 

XX ( IA) *XP 
YY (IA) *YP 

WRITE ( 6 , * )  'XP: • ,XP, 'YP: ' ,YP 
DADA=ADA/4. 

ADA1-ADA1+DADA 
30  CONTINUE 
C 

C******  GENERATING  POINTS  ON  RADIAL  PORTION  OF  RIM 
C 

DELTA-ASIN (XP/RR) 

RIM-RR-RTH 
DRTH-RTH/4 . 
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DX-RR+SIN ( DELTA) -RIM+SIN ( DELTA) 

DDX-DX/4 . 0 
DO  40  1-1,4 
IA— 24+ (1-1) 

XPRF ( I ) — XP-DDX 
YPRF ( I ) -YP-DRTH 
XX (IA) -XPRF (I) 

YY ( IA) —YPRF ( I ) 

DRTH-DRTH+RTH/4 . 

DDX-DDX+DX/4 . 0 

WRITE ( 6 ,  * )  ‘XP: ' , XPRF (I) , 'YP: • ,YPRF(I) 
40  CONTINUE 

XP— XPRF ( 4 ) 

YP— YPRF ( 4 ) 

C 

C******  GENERATING  POINTS  ON  INNER  RIM  SURFACE 
C 

DELTA-ATAN ( XP/ YP) 

XLEG- ( (XP*XP) + ( YP*YP) ) ** . 5 
DELTA1— DELTA/ 4 . 0 
DO  50  1-1,4 
IA— 28+ (1-1) 

DELTA— DELTA-DELTA1 
XPRF ( I ) — XLEG+SIN ( DELTA) 

YPRF ( I ) —XLEG* COS ( DELTA) 

XX (IA) -XPRF (I) 

YY ( IA) —YPRF ( I ) 

WRITE ( 6 , * )  'XP: ' ,XPRF(I) , 'YP: ' ,YPRF(I) 
50  CONTINUE 

C 

C  ***************************************** 

C  *  THIS  SECTION  USES  COORDINATES  FROM  * 

C  *  THE  PREVIOUS  SECTIONS  TO  CALCULATE  * 

C  *  NODAL  NUMBERS  AND  COORDINATES .  THERE  * 

C  *  ARE  319  NODES  USED  IN  THIS  MODEL.  * 

C  ***************************************** 

c 

c******  GENERATING  NODES  1-121 
C 

II-O 

DO  60  1-1,11 
IA— 1+(I-1) *11 
YA-YY(I) 

XA-XX(I) 

AZ-XA* . 15 
BZ-XA* . 30 
CZ-XA* . 50 
DZ-XA* .75 
XN(IA) — XA 
XN ( IA+10) =XA 
XN ( IA+1 ) — XA+AZ 
XN ( IA+9 ) — XA-AZ 
XN ( IA+2 ) — XA+BZ 
XN ( IA+8 ) -XA-BZ 
XN(IA+3 ) — XA+CZ 
XN ( IA+7 ) — XA-CZ 
XN ( IA+4 ) — XA+DZ 
XN ( IA+6 ) — XA-DZ 
XN(IA+5)«0.0 
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60 

C 

C* 

C 


75 


70 


71 

C 


80 

C 


DO  60  J-1,11 
II-II+l 
YN(II)-YA 
CONTINUE 

*  GENERATING  NODES  122-209 

DANG-2 . 5 
SANG— 10.0 
DO  70  1—1,8 
DXANG— SANG 
DO  75  J=1 , 4 
XANG ( I ,  J ) —DXANG 
DXANG-DXANG-DANG 
WRITE(6, *)  XANG(I,J) 

XANG (I , J) —XANG (I , J) *PI/180. 

CONTINUE 
DANG— DANG+2 . 5 
SANG— SANG+10 . 

CONTINUE 
DO  71  1-1,8 
XANG (1,4) —0 . 0 
CONTINUE 

DO  80  1-1,8 

IA— 122+ (1-1) *11 
XA— XX(11+I) 

YA=YY(11+I) 

AZ— XA*. 15 
BZ— XA* .30 
CZ— XA* . 50 
DZ— XA* .75 
DZ— DZ* ( 1 . + (I* . 02) ) 

XN(IA) — XA 
XN ( I A+ 1 0 ) — XA 

XN ( IA+ 1 ) — XA+AZ *COS ( XANG ( I , 1 ) ) 

XN(IA+9) — XN (IA+1) 

XN (IA+2) — XN (IA+1) + (BZ-AZ) *COS (XANG (I , 2) ) 
XN ( IA+8 ) — XN ( IA+2 ) 

XN(IA+3) — XN (IA+2 ) + (CZ-BZ) *COS (XANG (I , 3) ) 
XN(IA+7) — XN(IA+3) 

XN ( IA+4 ) -XN ( IA+3 ) + ( DZ-CZ ) *COS (XANG ( I , 4 ) ) 
XN (IA+6) — XN(IA+4) 

XN(IA+5) — 0 . 0 
IF(I.EQ.7)  CZ-CZ*1.1 
IF(I.EQ.8)  CZ— CZ*1 . 15 
YN(IA) -YA 
YN(IA+10) — YA 

YN(IA+1) -YA-AZ*SIN(XANG (I , 1) ) 

YN(IA+9) -YN (IA+1) 

YN(IA+2) =YN(IA+1) - (BZ-AZ) *SIN (XANG (I , 2) ) 
YN ( IA+8 )-YN( IA+2) 

YN ( IA+3 ) -YN ( IA+2 )- (CZ-BZ )*SIN (XANG (I, 3) ) 
YN ( IA+7 ) — YN ( IA+3 ) 

YN ( IA+4 ) -YN ( IA+3 ) - ( DZ-CZ ) *SIN (XANG (1,4)) 
YN ( IA+6 )-YN( IA+4) 

YN(IA+5) — YN (IA+4) 

CONTINUE 
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DZ-DZ/1.16 

CZ-CZ/1.15 

C 

C******  GENERATING  NODES  210-233 
C 

DO  85  1-1,3 
J»20+(I-1) 

XINC(I) -ATAN(XX(J)/YY (J) ) 

85  CONTINUE 

DINC-XINC(l)/2. 

DO  90  1-1,3 

IA-210+(I-1)*8 

ya*yy(19+i) 

XA«XX(19+I) 

XN(IA)— XA 
XN(IA+7)-XA 

XN ( IA+1) -XN ( IA) +AZ*SIN ( XINC (I) ) 
XN(IA+6)— XN(IA+1) 

XN(IA+2)— XN(IA+1)+(BZ-AZ)  *SIN(XINC(I)  ) 
XN ( IA+5) — XN ( IA+2 ) 

XN ( IA+3 ) -XN ( IA+2 )  +  ( CZ-BZ ) *  SIN ( XINC (I) ) 
XN(IA+4)— XN(IA+3) 

YN(IA)=YA 

YN(IA+7)«YA 

YN ( IA+1 ) -YA-AZ  *COS ( XINC (I) ) 

YN ( IA+6 ) =YN ( IA+1 ) 

YN (IA+2 )-YN (IA+1) -(BZ-AZ)*COS (XINC (I)) 
YN ( IA+5 )«YN( IA+2) 

YN ( IA+3 ) =YN ( IA+2 ) - (CZ-BZ) *COS (XINC (I ) ) 
YN  ( IA+4  )  ®*YN  ( IA+3 ) 

90  CONTINUE 
C 

C******  GENERATING  NODES  234-241 
C 

DELTA=ATAN (XX ( 2  3 ) / YY (23)) 
ALPHA*DELTA-ATAN (XN (206) /YN (206) ) 

XN (241) —XX (23) 

YN(241)—YY (23) 

XN(234)=- XN(241) 

YN(234) — YN(241) 

XN(240) — XN(241) -AZ*SIN (DELTA) 

YN(240) =YN (241) -AZ*COS (DELTA) 

XN ( 2  3  5 ) =-XN (240) 

YN(235) — YN(240) 

XN(239) =XN(241) -BZ*SIN (DELTA) 

YN (239) -YN (241) -BZ*COS (DELTA) 

XN  (236)  — XN  (239) 

YN(236) — YN(239) 

XN (238) -XN (241) -CZ*SIN (DELTA) 

YN(238)— YN(241) -CZ*C0S (DELTA) 

XN(237) — XN(238) 

YN(237) -YN(238) 

C 

C******  GENERATING  NODES  242-319 
C 

RIM- ( YN ( 2  3  8 ) - YY (27)) /COS ( DELTA) 
DRIM-RIM/6 . 0 
DO  100  1-1,5 
IA— 238-(I-l)  *8 
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XADA ( I ) -ATAN (XN ( IA) / YN ( IA) ) 

XADA ( 6 ) * ATAN ( XN (205) /YN (205) ) 

100  CONTINUE 
C 

DO  110  1=1,6 
IA»238-(I-1)*8 
IX=254- (1-1) 

IF(I.EQ. 6)  IA=205 

XN ( 2 54 ) =XN ( 2 3 8 ) -DRIM*SIN ( XADA ( 1 ) ) 

YN ( 2  54 ) =YN (238) -DRIM*COS ( XADA ( 1 ) ) 
XRAD=XN (254) /SIN (XADA (1) ) 

IF(I.EQ.l)  GO  TO  111 
XN ( IX) =XRAD*SIN ( XADA ( I ) ) 

YN(IX) =XRAD*COS (XADA (I) ) 

111  XN (IX- ( 14-2*1) ) “-XN (IX) 

YN(IX-(14-2*I) )=YN(IX) 

110  CONTINUE 
C 

DO  119  J=l, 6 
DO  120  1=1,6 
IA=255+ ( J— 1)  *13+ (I— 1) 

XN ( IA) =XN ( IA-13 ) +DRIM*SIN (XADA (I) ) 

YN ( IA) =YN ( IA-13 ) -DRIM*COS (XADA (I) ) 
XN(IA+ ( 14-2*1) ) =-XN (IA) 

YN (IA+ ( 14-  2*1) )=YN(IA) 

120  CONTINUE 

119  CONTINUE 

DO  125  1=1,6 

IA=248+( (1-1) *13) 

XN(IA) =0 . 0 
YN(IA) =YN(IA-1) 

125  CONTINUE 
C 

DO  130  1=1,319 

WRITE (11, 1000)  I , XN ( I ) , YN ( I ) 

130  CONTINUE 

1000  FORMAT ( 110, 7X, • 0. 0 • , F10. 5, F10. 5) 

C 

C******  THIS  ENDS  THE  NODE  GENERATION  SECTION 

C 

C 

C  A**************************************** 

C  *  THIS  SECTION  GENERATES  THE  ELEMENTS  * 

C  *  AND  ELEMENT  NUMBERS.  TWO-DIMENSIONAL* 

C  *  4 -NODE D  ELEMENTS  ARE  USED.  FOR  THIS  * 

C  *  MODEL  THERE  ARE  276  ELEMENTS.  * 

C  ***************************************** 

c 

c******  GENERATING  ELEMENTS  1-183 
C 

11=0 

DO  140  1=1,18 
IA=12+(I-1) *11 
IB=13+ (1-1) *11 
IC=2+ (1-1) *11 
ID=1+ (1-1) *11 
DO  145  J-1,10 
II-II+l 
LA=IA+ (J-l) 


LOIC+(J-l) 

LD=ID+ ( J-l) 

145  WRITE(11, 112)  II, LA , LB , LC , LD 

140  CONTINUE 

WRITE (11, 112) 181,210,211,200,199 
WRITE (11 #112 >182,211, 212, 201, 200 
WRITE (11, 112) 183,212,213,202,201 
C 

C******  GENERATING  ELEMENTS  184-195 
C 

11=183 

DO  150  1=1,3 
IA=218+ (1-1) *8 
IB=219+ (1-1) *8 
IC=211+ (1-1) *8 
ID=210+(I-1) *8 
DO  155  J=1 , 3 
11=11+1 
LA=IA+ (J-l) 

LB=IB+ (J-l) 

LC=IC+ (J-l) 

LD=ID+ (J-l) 

155  WRITE (11,112)11, LA, LB, LC, LD 

150  CONTINUE 

WRITE (11, 112) 193,214,215,207,206 
WRITE (11, 112) 194,215,216,208,207 
WRITE (11, 112) 195,216,217,209,208 
C 

C******  GENERATING  ELEMENTS  196-203 
C 

11=195 

DO  160  1=1,3 
IA=222+ (1-1) *8 
IB=223+ ( I — 1 ) *8 
IC=215+ (1-1) *8 
ID=214+ (1-1) *8 
DO  165  J=l, 3 
11=11+1 
LA=IA+ (J-l) 

LB=IB+ (J-l) 

LC=IC+ (J-l) 

LD=ID+ (J-l) 

165  WRITE(11, 112)11, LA, LB, IC,LD 
160  CONTINUE 
C 

C******  GENERATING  ELEMENTS  204-213 
C 

11=204 

DO  170  1=1,3 
LA=242+ ( 1-1) 

LB=243+ (1-1) 

LC=229- ( 1-1) *8 
LD=237- (1-1) *8 
11=11+1 

170  WRITE(11, 112)11, LA, LB, LC,LD 

WRITE (11, 112) 208, 245, 246, 202, 213 
WRITE (11, 112) 209, 246, 247, 203, 202 
WRITE (11, 112) 210,247,248,204,203 


m\j.  it.  205,  204 

WRITE (11, 112 >212,249,250 ,206, 205 

WRITE (11, 112) 213, 250, 251, 214 ,206 
C 

C******  GENERATING  ELEMENTS  214-276 
C 

11=213 

DO  180  1=1,3 
LA=251+ (1-1) 

LB=252+ (1-1) 

LC=222+ (1-1) *8 
LD=214+(I-1) *8 
11=11+1 

180  WRITE(11, 112)11, LA, LB, LC,LD 

11=216 

DO  190  1=1,5 

IA=255+ (1-1) *13 
IB=256+ (I-l)*13 
IC=243+ (I-l) *13 
ID=S242+  (1-1)  *13 
DO  195  J-l, 12 
11=11+1 
LA=IA+(J-1) 

LB=IB+ (J-l) 

LC=IC+ ( J-l) 

LD=ID+ (J-l) 

195  WRITE(11, 112)11, LA, LB, LC,LD 

190  CONTINUE 
C 

112  FORMAT (15 , 4X, '8' ,4I5,23X, '21' ,4X, '1*) 

C 

C******  THIS  ENDS  THE  ELEMENT  GENERATION  SECTION 
C 

STOP 

DEBUG  SUBCHK 


V 


oooooooooooftftnoftoftnnftoooooooooooooo 


C***************************** *************** ****** 


* 

VARIABLE  SPEED/LOAD  GENERATION  PROGRAM  * 

* 

THIS  PROGRAM  GENERATES  THE  SAP6  TIME  * 

FUNCTIONS  DESCRIBING  THE  LOAD  MOVING  * 

OVER  THE  INVOLUTE  OF  A  GEAR  TOOTH.  THE  * 

SPEED  OF  THE  MOVING  LOAD  ON  THE  INVOLUTE  * 

OF  A  TOOTH  VARIES  AS  * 

* 

V(T) — RB*OMEGA**2*T  * 

* 

WHERE  RB  IS  THE  BASE  RADIUS  OF  THE  GEAR  * 

OMEGA  IS  THE  ANGULAR  VELOCITY  OF  THE  GEAR  * 

AND  T  IS  THE  ELLAPSED  TIME.  THE  LOAD  * 

FUNCTIONS  WERE  TAKEN  FROM  WALLACE ,  AND  ARE  * 

* 

F (T) =F0 ( 1-COS (T/ALPHAT) ) /2  0<T<ALPHAT  * 

* 

F(T)-FO  ALPHAT<T<( 1-ALPHA) TF  * 

* 

F (T) -F0 ( 1-COS (T-ALPHATF) /ALPHATF) )/2  * 

( 1 -ALPHA ) TF<T<TF  * 

* 

WHERE  TF  IS  THE  TOTAL  TIME,  AND  ALPHA  IS  * 

A  FACTOR  DEPENDENT  ON  THE  CONTACT  RATIO.  * 

* 

THE  USER  MUST  ENTER  THE  FOLLOWING  PARAMETERS  * 

* 

PRESSURE  ANGLE-PANG  * 

PITCH  RADIUS  -RP  * 

ADDENDEM  =AD  * 

CIRCULAR  PITCH=CIRP  * 

BACKLASH  -BACKL  * 

MAX  VELOCITY  -VMAX  * 


* 

C************* ******************************* ****** 

DIMENSION  XX (40) ,YY(40) , FORCE (11 , 500) ,TIME(500) ,T(500) 

DIMENSION  XINC(ll) ,DXINC(11) ,FX(11,500) ,FY(11,500) ,XN(15) ,YN(15) 

PI-3.141592654 

READ (  5 ,  * )  PANG 

PANG— PANG* PI/ 180 . 

READ(5, *)  RP 
READ(5, *)  AD 
READ (  5 ,  * )  CIRP 
READ ( 5 ,  * )  BACKL 
READ(5, *)  VMAX 
RB— RP*COS (PANG) 

THP— ( (RP*RP-RB*RB) / (RB*RB) ) **0. 5 
CIRTH— CIRP/2 . -BACKL 
AL-THP-ATAN(THP) +CIRTH/ (RP*2 . ) 

RO-RP+AD 

C 

C ************************************************** 

C  THE  INVOLUTE  COORDINATES  AND  INVOLUTE  NORMALS  * 

C  ARE  CALCULATED  IN  THIS  SECTION.  * 

C ************************************************** 

c 

DO  5  1-1,11 
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O  Ul  to  H*  *  *  *  * 


IA=1+(I-1) 

HH=(RO-RB)/10. 

H»R0-HH*(I-1) 

ANG=( (H*H-RB*RB) / (RB*RB) ) **0. 5 
Y=RB* (COS (ANG) +ANG*SIN (ANG) -1 . ) 

X— RB*  (SIN  (ANG)  -ANG*COS  (ANG)  ) 

XP=RB*SIN ( AL) +X*COS (AL) +Y*SIN (AL) 

YP=RB*COS (AL) -X*SIN (AL) +Y*COS (AL) 

XX(IA)=XP 

YY(IA)=YP 

XN (I) =-SIN (AL) *SIN (ANG) -COS (AL) *COS (ANG) 

YN(I) =SIN (AL) *COS (ANG) -COS ( AL) *SIN (ANG) 

CONTINUE 

************************************************* 

THE  DISTANCE  ALONG  THE  INVOLUTE  IS  DETERMINED  * 

FROM  WHICH  THE  TOTAL  TIME  AND  ANGULAR  VELOCITY  * 

ARE  FOUND.  * 

************************************************* 

S=0.0 

DO  6  1=1,10 

XL=XX(I+1) -XX (I) 

YL=YY(I)-YY(I+1) 

XINC(I)=(XL**2+YL**2) **0.5 
S=S+XINC(I) 

WRITE ( 6 , * )  ' XINC (I) =  ',XINC(i),'  S=  *,S 

CONTINUE 
TF=2 . *S/VMAX 

XOMEGA* ( VMAX/ (TF*RB) ) ** . 5 

WRITE ( 6 , * )  *  TF=  ' , TF , '  OMEGA=  • , XOMEGA 
TALPHA=0.28*TF 

WRITE ( 6 , * )  ' TALPHA™  ' ,TALPHA, ’  ( 1-ALPHA) *TF=  ' , (1. 28) *TF 

************************************************* 

THE  DISTANCE  BETWEEN  EACH  NODE  IS  DIVIDED  UP  * 

INTO  20  EQUAL  INCREMENTS  AND  THE  TIME  AND  * 

AND  FORCE  VALUES  ARE  DETERMINED  FOR  EACH.  * 

************************************************* 


DIST=0 . 0 
DO  10  1=1,10 

DXTNC ( I ) =XINC ( I ) / 2  0 . 

DO  15  J=l, 20 
JJ=J+20* (I-l) 

DIST-DIST+DXINC ( I ) 

TIME ( JJ) = ( 2 . *DIST/ (RB*X0MEGA**2 ) ) **0 . 5 
T(JJ)=TIME(JJ) 

IF (TIME (JJ) .GT.TALPHA)  GO  TO  11 

FORCE ( I , J) = (FO/2 . ) * ( 1 . -COS (PI*TIME ( JJ) /TALPHA) ) 

GO  TO  15 

IF(TIME(JJ) .GT. (TF-TALPHA) )  GO  TO  12 
FORCE ( I , J ) =FO 
GO  TO  15 

FORCE (I ,J) = (FO/2 . ) * ( 1 . -COS (PI* (TF-TIME (JJ) ) /TALPHA) ) 
CONTINUE 
CONTINUE 

DO  20  1=1,1 
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DO  30  J-1,20 
JJ-J+20*(I-1) 

FX ( I , J+20) -FORCE (I, J) * ( (XINC ( I ) - (RB* (XOMEGA**2 ) /2 . ) 
$*TIME(JJ)**2)/XINC(I) )*XN(I) 

FY (I /J+20) -FORCE (I ,J) * ( (XINC (I) - (RB* (XOMEGA**2)/2 . ) 

$*TIME(JJ) **2) /XINC(I) ) *YN(I) 

FACT=(RB* (XOMEGA**2)/2 . ) * (TIME(JJ+180) **2-TIME (180) **2) 
FX ( 11 , J) —FORCE ( 10 , J) * ( FACT/XINC (10) ) *XN(11) 

FY(11, J) -FORCE (10, J)*(FACT/XINC(10))*YN(11) 

30  CONTINUE 

20  CONTINUE 

DO  25  1-2,10 
DO  35  J— 1,20 
JJ— J+20* (1-2) 

XTIME1— TIME( (1-1) *20) 

IF(I.NE.2)  GO  TO  36 
XTIME-0 . 0 
GO  TO  37 

36  XTIME— TIME ( (1-2) *20) 

37  FACT=( (RB* (XOMEGA**2) )/2 . ) * (TIME(JJ) **2-XTIME**2) 

FACT1— ( (RB* (XOMEGA**2) )/2 . ) * (TIME(JJ+20) **2-XTIMEl**2) 

FX ( I , J ) -FORCE (I-1,J)*( FACT/XINC ( I- 1 ) ) *XN(I) 

FY ( I , J) -FORCE (1-1, J) * ( FACT/XINC ( 1-1) ) *YN (I ) 

FX(I , J+20) —FORCE (I ,J) * ( (XINC (I) -FACT1)/XINC(I) ) *XN(I) 

FY(I, J+20) —FORCE (I, J) * ( (XINC (I) -FACT1)/XINC(I) ) *YN(I) 

35  CONTINUE 

25  CONTINUE 
C 

C** ************************************** ********* 

C  THE  REMAINDER  OF  THE  PROGRAM  WRITES  THE  TIME  * 

C  FUNCTION  VALUES  (TIME, FORCE)  TO  A  FILE  FOR  * 

C  USE  IN  A  FINITE  ELEMENT  CODE.  THIS  PROGRAM  * 

C  IS  FORMATTED  FOR  SAP6.  * 

C************************************************* 

c 

NFN-22 

WRITE (11, 99)  NFN 

99  FORMAT (3X, 12, 4X, • 0 • , 4X, ' 0 ' , 3X, 'NT' ,4X, *1' ,6X, ' DPER' ,7X, *0.0') 
IFN-0 

DO  60  1-1,11 
NP— 11+11* (1-1) 

IC-2 

DO  60  J— 1 , 2 
IFN-IFN+1 

WRITE(11, 98)  NP, IC, IFN 
IC-3 

98  FORMAT (15 , 4X, Il,I5,4X,'l', 7X, '1.0' ,24X, '0') 

60  CONTINUE 

WRITE(11,97) 

97  FORMAT (/•  •) 

C 

C*****  TIME  FUNTION  DATA  **** 

C 

TF-1.0 

TO-O.O 

FO-O.O 

C********  TIME  FUNCTION  FOR  NODE  11  ******* 

C 


NLP-24 
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WRITE (11 ,101)  NLP 

101  FORMAT (3X, 12 , 11X, 'TIME  FUNCTION' , 50X, • 1 ') 

100  FORMAT (3 (E12 . 6, F12 . 0) ) 

WRITE ( 11 , 100)  T0,F0,T(1)  ,FX(1,21)  ,T(2)  ,FX(1,22) 

DO  200  1-3,18,3 

200  WRITE (11, 100)  T(I) ,FX(l,I+20) ,T(I+1) , FX(1, 1+21) ,T(I+2) , FX(1, 1+22) 
WRITE (11, 100)  TF, F0,TF, F0,TF,F0 

WRITE (11, 101)  NLP 

WRITE (11, 100)  TO ,F0,T(1) , FY (1,21) ,T(2) ,FY(1,22) 

DO  201  1=3,18,3 

201  WRITE(11, 100)  T(I) ,FY (1,1+20) ,T (1+1) , FY ( 1 , 1+21) ,T(I+2) ,FY(l,I+22) 
WRITE (11, 100)  TF , FO , TF , FO , TF , FO 

C******  TIME  FUNCTIONS  FOR  NODES  22  THRU  110  ***** 

C 

DO  202  1=2,10 
NLP=45 

WRITE (11, 101)  NLP 
K=20*(I-2) 

T(0)=0.0 

WRITE (11, 100)  T0,F0,T(K) , F0,T(1+K) , FX(I, 1) 

DO  203  J=2 ,38,3 

203  WRITE (11, 100)  T(J+K) ,FX(I,J) ,T(J+1+K) , FX(I, J+l) ,T(J+2+K) ,FX(I, J+21 
WRITE (11, 100)  TF , FO , TF , FO , TF , FO 
WRITE (11, 101)  NLP 

WRITE (11, 100)  TO, F0,T(K) ,FO,T(l+K) ,FY(I,1) 

DO  205  J=2 ,38,3 

205  WRITE (11 , 100)  T(J+K) ,FY(I,J) ,T(J+1+K) ,FY(I, J+l) ,T(J+2+K) ,FY(I,J+2; 
WRITE (11, 100)  TF , FO , TF , FO , TF , FO 

202  CONTINUE 
C 

C****  TIME  FUNCTION  FOR  NODE  121  ******* 

C 

NLP=24 

WRITE (11, 101)  NLP 

WRITE (11 , 100)  TO,FO,T(180) ,F0,T(181) ,FX(11,1) 

DO  206  1=2,17,3 
K—I+180 

206  WRITE (11, 100)  T(K) , FX(11, I) ,T(K+1) , FX(11, 1+1) ,T(K+2) ,FX(ll,I+2) 
WRITE (11, 100)  T(200) ,FX(11,20) ,TF,F0,TF,F0 

WRITE (11, 101)  NLP 

WRITE (11, 100)  T0,F0,T(180) ,F0,T(181) ,FY(11,1) 

DO  207  1=2,17,3 
K—I+180 

207  WRITE (11, 100)  T(K) , FY ( 11 , I) ,T(K+1) ,FY(11,I+1) ,T(K+2) ,FY(ll,I+2) 
WRITE (11, 100)  T(200) ,FY(11,20) ,TF,F0,TF,F0 

C 

C****  ECHO  CHECK  ****** 

C 

DO  40  1=1,11 

IF(I.NE.l)  GO  TO  52 
DO  42  J—l, 20 
JJ=J+20*(I-1) 

WRITE (12 ,51)  T(JJ) ,FX(I,J+20) 

42  CONTINUE 

GO  TO  40 

IF(I.EQ.ll)  GO  TO  53 
DO  45  J—l, 20 
JJ— J+20* (1-2) 

WRITE(12 , 51)  T(JJ) ,FX(I,J) 
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CONTINUE 
DO  46  J-1,20 
JJ-J+20*(I-1) 

WRITE (12 , 51)  T(JJ) ,FX(I,J+20) 
CONTINUE 
GO  TO  40 
DO  50  J“l,20 
JJ=J+20*(I-2) 

WRITE (12 , 51)  T ( JJ) / FX (I , J) 
CONTINUE 
CONTINUE 

FORMAT ( F12 . 8 , 12X , F8 . 0 ) 

S' 


APPENDIX  2 

Contact  Point  Velocity 
Of  Meshing  Spur  Gear  Teeth 


Por  a  meshing  spur  gear  pair,  the  speed  of  the  contact 
point  varies  with  time  as  it  moves  from  the  tip  of  the  tooth  to 
the  base  circle. 

From  equations  (2-3)  and  (2-4)  and  Figure  A2-1  the  coor¬ 
dinates  of  a  point  Bi  on  the  involute  are  given  by; 

XBj.  -  -RB(sin  0i  -  6iCOS6i)  (A2-1) 

YBi  ■  RB(COS01  ♦  QisinQi  -  1)  (A2-2) 


Differentiating  (A2-1)  and  (A2-2)  with  respect  to  0j_; 

dXB  -  __  „  .  „ 

1=  -  RB  0j  sm0i 

dYBi 

"  de:"  ~  RB0iCOS0i 
i 

and  combining  to  give  the  resultant; 


,di=,2  dXBi,2  .  ,dYBi.2 


RB2  0?  (sin20 j+cos20i) 


(A2-3) 

(A2-4) 


dr 


d0. 


RB0  j 


(A2-5) 


Multiplying  through  by  d  0i  ■  §i ,  realizing  that  0i  and  r  vary 
with  time,  we  have 

at  -  V  -  RB0i0i 

And  if  we  define  0i  ■  wt; 


V(t)  -  RBw2t 


(A2-6) 


where; 


\ 


Figure  A2-1:  Contact  point  goemetry 
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V(t)  *  the  speed  of  the  contact 
point  Bj,  at  time  t 


w  *  angular  velocity  of  the  gear  (rad/sec) 


At  t  *  0.0,  corresponding  to  the  load  at  the  tip  of  the  tooth, 
the  speed  is  zero,  and  at  t  »  TF,  time  corresponding  to  the 
load  positioned  at  the  base  circle,  the  speed  is  maximum. 

The  position  of  the  load  as  a  function  of  time  is  found  by 
integrating  equation  (A2-6)  with  respect  to  time.  By  defining 
the  speed  as  the  first  derivative  of  the  position; 


V(t) 

dS(  t) 

S(t) 


If  (t>  *  RBu^t 
at 


RBoj2  J  t  dt 
^RBu?t2  +  c 


(A2-7) 


Evaluating  the  constant  of  integration  for  t  »  Tj.,  the  time  at 
the  tip  of  the  tooth; 


then; 


S(ti)  -  0.0  -  %RBu>2Ti2+  C 


C  »  -3£RBu)2t| 

The  displacement  then  becomes; 


S(t)  -  3£RBo2t2  -  ^RBa)2Ti2 


(A2-8) 


(A2-9 ) 


If  we  assume  that  the  time  at  the  tip  position  is  zero  the 
constant  term  vanishes  and  we  are  left  with  only  the  first  term 


on  the  right  hand  side  of  equation  A2-9  defining  the  position 
along  the  involute. 
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APPENDIX  3 

Deformed  Shapes  of  Gear  Tooth 

1.  Static  Loading 

2.  Modal  Analysis 

3 .  Dynamic  Response 


STATIC  DEFLECTION  ANALYSIS:  GEAR  #1  (BEAM  CONSTR) 
UNDEFORMED  SHAPE 


Figure  A3-1 
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STRTIC  DEFLECTION  RNRLYSIS :  6ERR  *1  (BERM  CONSTR) 
STRTIC  LORD  CRSE  1  t  =  0.1 


SEPTEMBER  09.  1983  01:37:07 
1  RXIS”3  RLPHfi=  0.00  BETR=  0.00 

DEFLECTION  SCALE  FflCTCRMH  . 1 


Figure  A3-2 


STATIC  DEFLECTION  ANALYSIS:  GEAR  *1  (BEAM  CONSTR) 
STATIC  LOAD  CASE  2  t  =  0.3 

SEPTEMBER  09.  1983  01  *37  *07  ] 

IflXl S=»3  ALPHA-  0.00  BETA-  0.00 
DEFLECTION  SCALE  FACTCKW8.6 

> 

> 


Figure  A3-3 


STATIC  DEFLECTION  ANALYSIS:  GEAR  «1  (BEAM  CONSTR) 
STATIC  LOAO  CASE  3  t  -  0.5 


9EPTEMKA  88.  19(3  BlUTigr 

IflKlS-3  BLPHfV  0.88  KTA-  0.08 

deflection  scale  factomt.i 


Pigurt  A3-4 


STATIC  DEFLECTION  ANALYSIS:  GEAR  #1  (BEAM  CONSTR) 
STATIC  LOAD  CASE  4  T  =  0.8 

SEPTEMBER  OS.  1983  0  H37I07 

IRXIS-3  ALPHA-  0.00  BETA-  0.00 

DEFLECTION  SCALE  FBCJTWH  . 


Figure  A- 5 


STRTIC  DEFLECTION  RNALYSIS:  6ERR  *1  (RIM  INCLUDED) 
STATIC  LORD  CASE  1  t  =  O.l 

SEPTEMBER  10-  1983  00«10i05  7 

IRXIS-3  RLPHR=  0.00  BETR*  0.00 

DEFLECT I  ON  SC RLE  FRCTORW3.8 

X 


Figure  A3-6 


STRTIC  DEFLECTION  RNRLYSIS:  GERR  #1  (RIM  INCLUDED) 
STRTIC  LORD  CRSE  2  t  =  0.3 


SEPTEMBER  10.  1383  00  *  10  >05 

1AXIS-3  ALPHA-  0.00  BETA-  0.00 

DEFLECTION  SCALE  FACTCRM1  .6 


Figure  A3-7 
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STATIC  DEFLECTION  ANALYSIS 
STATIC  LOAD  CASE  3  t  =  0.5 

SEPTEMBER  10.  1983  00  1  10  105 

IRXIS-3  RLPHR-  0.00  BETA-  0.00 

DEFLECTION  SCALE  FACTCJS67  .S 


GEAR  #1  (RIM  INCLUDED) 


■ 

111 


BhbB 


Figure  A3 -8 


STATIC  DEFLECTION  ANALYSIS:  GEAR  #1 
STATIC  LOAD  CASE  4  t  =  0.8 

EPTEMBER  10.  1983  OOllOlOS 

AX IS- 3  ALPHA-  0.00  BETA-  0.00 

CFLECTION  SCALE  FACTCf&C*  .3 


(RIM  INCLUDED) 


Figure  A3-9 


Modal  Analysis 
Timoshenko  Beam  Constraints 
Rim  Constraints 
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MODAL  ANALYSIS,  PLANE  STRAIN  ELEMENT:  GEAR  #1 
UNDEFORMED  SHAPE 


Ut«  1G.  1983  00  : 10  i  15* 

IflXl'3=3  RLPHfl=  0.00  BETfl=  0.00 


Figure  A3-1Q 


MODAL  ANALYSIS,  PLANE  STRAIN  ELEMENT:  GEAR  #1 
MODE  1  Freq.=  ,9101  E5  (CPS) 


UNE  10.  1983  19  *32  : 3 9 * 

inxis-3  ALPHA-  0.00  BETA”  0.00 

DEFLECTION  SCALE  FACTOR”  0  .21278C 


MODfiL  ANALYSIS ,  PLANE  STRAIN  ELEMENT:  GEAR  #1 
MODE  2  Freq.=  .2159  E6  (CPS) 


Figure  A3-12 
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MODRL  RNRLYSIS,  PLRNE  STRRIN  ELEMENT:  GERR  #1 
MODE  3  Freq.=  .2289  E6  (CPS) 

UNE  10.  1083  19  *32  :39' 

IPXISO  OLPHR=  O.On  BEl'H-  n  .00 
DEFLECTION  SCALE  FACTOR=  0  .21756-E 


Figure  A3-13 
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MODAL  ANALYSIS.  PLANE  STRAIN  ELEMENT:  GEAR  # 1 
MODE  4  Freq.=  .4218  E6  (CPS) 

UNE  10.  1983  19  :32  :39  ’ 

I  flX  I  S=  3  RL  PHP=  0.00  BET  0=  0.00 

DEFLECTION  SCFiLE  FflCTOR=  0  . 2  1 18 *£ 


Figure  A3-14 


MODAL  ANALYSIS,  PLANE  STRAIN  ELEMENT:  GEAR  #1 
MODE  5  Freq.=  .4794  E6  (CPS) 

'JNC  10.  1383  19  1 3 2  1 3 9  * 

! AX ! C " 3  RLAHfl=  n,00  BETR=  0.00 

DEFLECTION  SCALE  FRCTOR=  0  .29789-E 


Figure  A3-15 


MODfiL  ANALYSIS.  PLANE  STRAIN  ELEMENT:  GEAR  #1 
MODE  1  Freq .  =  .6545  E5  (CPS) 

UNE  16.  1983  00  i  10  s  15 *  7 

I  AXIS’*  3  ALPHA*  0.00  BETA*  0.00  ! 

DEFLECTION  SCALE  FACTO?*  0  .27565E  j 

* - Y 


Figure  A3-16 
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MODAL  ANAlySIL  ->lrtNE  S T RP [ N  ELEMENT:  GEAR  «1 
MODE  3  Freq.-  .1517  E6  (CPS) 


JNF  M  l«»f  0 C  :  ■  0  !  S ' 

' »>  i  S  «  i  Rt  ««.  5  .  f  -J  yf  r  f  i,  0  .  3  U 

OCFi.ttTlOfc  SC  Hit  iHtrj).  u.ih'H! 


I 


? 

I 

X-  -Y 


Figure  A3-18 
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MODRL  RNRLYSIS.  PLRNE  STRRIN  ELEMENT 
MODE  4  Freq.=  .2034  E6  (CPS) 


GERR  #1 


UNE  16  1983  00  s  10  !  15 * 

IAXIS-3  ALPHA-  0.00  &ETA- 
OEFlECTI  ON  SCALE  FACTOR-  0  .4  2  67  7-E 


mm 

■Ml 


Figure  A3-19 


MODAL  ANALYSIS,  PLANE  STRAIN  ELEMENT 
MODE  5  Freq .  =  .2319  E6  (CPS) 

UNE  16.  1  983  00  :  10  i  15 ' 

I  AXIS”  3  ALPHA-  0.00  BETA-  0.00 

OEFLECTION  SCALE  FACTOR-  0  .5253  8-E 


Dynamic  Response 


Gear  Tooth  Profile  Deflections 


In  the  nodal  response  analysis,  nodal  superposition  is 
used  to  sun  the  effects  of  all  the  included  nods  shapes  pro¬ 
ducing  the  desired  response.  The  nore  nodes  included,  the  sore 
accurate  the  analysis.  However,  there  is  a  practical  linit  to 
the  nunber  of  nodes  used  dictated  by  the  type  of  loading  used 
(inpact,  constant,  sinusoidal,  etc.)  and  of  course  con- 
putational  efficiency.  As  previously  nentioned,  the  first  ten 
siodes  are  used  for  this  analysis. 

To  better  conprehend  the  dynanic  deflection  phenosMoa,  the 
deflected  state  of  the  tooth  profile,  subject  to  the  inpact 
loading  case,  is  plotted  for  ditferent  load  positions  with 
v*«0.01.  Figure  A3-21  shows  the  under forned  shape  and  Figures 
A3-22  through  A3-25  are  the  various  deforned  shapes.  Froa  the 
figure(s),  one  can  see  those  nodes  which  contribute  noticeably 
to  the  overall  deflection  of  the  tooth.  It  is  apparent  that 
only  the  first  three  or  four  have  a  najor  effect  (see  Figure 
A3- 19,  T«0 . 5 ) .  Thus  the  problen  of  local  defornation  encoun¬ 
tered  in  the  static  analysis  is  not  apparent  here. 
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DYNAMIC  DEFLECTIONS  OF  GEAR  TOOTH  WITH  MOVING  LOADS 


Figure  A3-22 


DYNAMIC  DEFLECTIONS  OF  GEAR  TOOTH  WITH  MOVING  LOADS 


I 

! 


Figure  A3-23 


V*=0.01  T=0.5 

DYNAMIC  DEFLECTIONS  OF  GEAR  TOOTH  WITH  MOVING  LOADS 


Figure  A3-24 
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.  V  '  ’  > 


V*=0.D1  T=0.8 

DYNAMIC  DEFLECTIONS  OF  GEAR  TOOTH  WITH  MOVING  LOADS 


Figure  A3-25 


WS  PLACE  MENT(INCHES)  DtSPLACEMENT(INCHES)  DISPLACEMENT  (INCHES) 


«*OStiiON 


VEL=  10.0  IN/SEC  .  MASS  =  1.0  .  BF AM  LENi.i‘<--  w  *,  N 


DYNAMIC  RESPONSE  OF  MESHING  CANTILEVER  BEAMS,  WITH  INf  WTia  = 


Figure  A4-16 
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APPENDIX  5 

Dm  of  Natural  Nodes  in  Equation 
of  Notion  for  Cantilever  B«tu 

1)  Constant  of  Integration 

2)  Evaluation  of  I2FBI1  AMD  I2PHI2 

3)  Derivatives  of  the  Node  Shape  With  Time 
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1.)  Constant  of  Integration 

In  order  to  simplify  the  solution  of  the  equation  of 
motion  (5-26),  the  constant  of  integration  associated  with  the 
mode  shape  of  a  cantilever  beam  is  determined  such  that; 


p(  <TU)d£  *  1  (A5-1) 

Jo 

Since  the  beams  are  identical  and  the  integral  is  over  the 
length  of  the  beam,  the  constant  is  the  same  for  both  beams. 

The  first  natural  mode  of  vibration  for  a  uniform  fixed- 
free  cantilever  beam  is; 


where ; 


<t>(5)  ■  Ct  (sin 3L-sinh3L Main 8£-sinhB£) 


+  (cosgL+coshBL)  (cosB€-cosh0€) 


w2f 


Equation  (A5-2)  is  simplified  by  rewriting  with; 


SI  -  sin  BL-sinh  6L 


S2  «  cos3L+cosh$L 


We  then  have; 


<t>(  C>  ■  C[Slsinx-Slsinhx+S2cosx-S2coshx] 


where  x  is  used  in  place  of  the  argument  . 


(A5-2) 


(A5-3) 


209 


Squaring  equation  (A5-3); 

$2(5)  -  C2[Sl2(sin2x+sinh2x-2sinxsinhx) 

+  S22(cos2x+cosh2x-2cosxsinhx)  (A5-4) 

+2SlS2(8inx  cosx-sinx  coshx-cosx  sinhx+sinhx  coshx)] 

We  can  then  solve  for  the  constant  C  with  the  following 
substitution; 


C  =(____! - )**  (A5-5) 

p/  G(£)d£ 

0 

where  <j>  (  K  )*CG(  £ ) .  Whenever  the  relation  represented  by 
equation  (A5-1)  occurs  in  equations  (5-26)  it  is  replaced  by  1. 
Throughout  the  remainder  of  equations  of  motion  the  constant  of 
integration  C  is  represented  by  equation  (A5-5).  This  relation 
is  evaluted  using  Simpson's  Rule. 


2.)  Evaluation  of  I2PHI1  and  I2PHI2 

Using  the  arbitrary  constant  of  integration ,  C,  determined 
in  the  previous  section,  the  relationship; 

I2PHI1  -  I2PHI2  «f  $2(£)d€  (A5-6) 

J  0 

can  be  evaluated.  The  integrand  represents  the  second  deriva¬ 
tive  of  the  mode  shape  with  respect  to  the  length  variable, £  , 
quantity  squared.  The  first  derivative  of  the  mode  shape  with 
x  again  used  in  place  of  BO 


||  -  CftSl 


cosx-S2coshx-S2sinx-S2sinhx] 
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and  the  second  derivative; 


£—1  -  C82[-Slsinx-Slsinhx-S2cosx-S2coshx] 
d£2 

Squaring  this  relation  and  simplifying  yields  the  desired 
result; 

(d2^k  2  _  c2g4 [Sl(sin2x+sinh2x+2sinx  sinhx) 

+S22(cos2x+cosh2x+2cosx  coshx)  (A5-7) 

+2SlS2(sinx  cosx+sinx  coshx+cosx  sinhx+sinhx  coshx)] 

This  relation  is  then  used  in  equation  (A5-6)  which  is  eva¬ 
luated  using  Simpson's  Rule. 

3.)  Derivatives  of  the  Mode  Shape  with  Time 

When  the  constraint  equation  is  differentiated,  first  and 
second  derivatives  of  the  mode  shape  result  where  the  dependent 
variable,  £  ,  must  be  considered  as  a  function  of  time.  The 
first  derivative  of  the  mode  shape  with  respect  to  time  is; 

d<b(£;) 

1  »  CB£ [Slcosx-Slcoshx-S2sinx-S2sinhx]  (A5-8) 

where  x  is  used  in  place  of  the  argument  (J£  and  £  is  the  velo¬ 
city  of  the  moving  load.  Taking  the  second  derivative,  we 
have; 

j2 , if i 

— *  CtSl(-B2  £2(sinx+sinhx)  4rV( cosx-coshx ) ) 
dr 

-  S2( 82£2<cosx+cosh) 

+3£ (sinx+sinhx) ) ]  (A5-9) 

where  £  is  the  acceleration  of  the  moving  load.  For  a  moving 
load  with  constant  speed,  £  is  of  course  zero. 


APPENDIX  6 


Development  and  Solution  of 
Equations  of  Motion  with 
Inertial  Terms  Removed 


By  eliminating  the  inertial  terms  from  the  first  four 
equations  of  (5-26)#  and  using  the  undifferentiated  form  of  the 
constraint  equation  (5-18)#  a  new  set  of  equations  for  the 
beam-mass  system  take  the  form; 


[A1(X}  -  {B} 


(A6-1) 


where ; 


[A] 


( Ml+MBl ) 
0 
0 
0 
0 


(M2+MB2 ) 
0 
0 
0 


0 

0 

I2PHI1 

0 

-PHI1 


0 

0 

0 

I2PHI2 

-PHI2 


1 

-1 

PHI1 
PHI  2 
0 


{X} 


(B} 


The  initial  conditions  for  this  system  are  determined  by 
arbitrarily  chosing  three  of  the  four  unknowns  composing  the 
constraint  equation.  Letting  the  beams  assume  static  deflec¬ 
tions  at  time  equal  to  zero#  and  XI  equal  to  zero;  the  initial 
conditions  are; 

P.PHI1 


qKO)  -  j2phi2 


q2  (0) 


P*PHI2 


I2PHI2 


xno)  -  o.o 


with; 


1 
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X2  ( 0 )  -  XI  +  qlPHIl  +  q2PHI2 


Solving  for  (X)  in  equation  (A6-1)  by  inversion; 

{X}  ■  [A]’^!} 

gives  the  values  of  the  unknowns  in  terms  of  XI  and  X2.  Using 
these  values  for  (X}»{X1  X2  ql  q2  X  }T,  the  system  is  integrated 
using  a  Runge-Kutta  integration  algorithm  solving  for  XI  and 
X2  and  X.  Then  equations  3  and  4  of  (A6-1)  are  used  to  solve  for  ql 
and  q2. 

Tests  were  done  for  constant  velocity  moving  loads  of  1.0r 
5.0,  10.0,  20.0  and  40.0  inches  per  second.  Plots  of  X2-X1, 
qlPHIl  and  q2PHI2  are  included  in  Figures  (A6-1)  through 
( A6-5 ) .  Comparing  these  to  Figures  (5-7)  to  (5-11)  show  that 
even  though  only  one  vibrational  mode  is  used,  the  results  are 
very  nearly  the  same  as  those  determined  for  the  massless 


beams 


DtSPLACEMENT(INCHES)  DISPLACEMENT (INCHES)  DISPLACEMENT(INCHES) 


TIME(SEC) 

VEL=  1.0  IN/SEC  ;  MASS=  1.0  ;  BEAM  LENGTH  =  0.5  IN 


DYNAMIC  RESPONSE  OF  MESHING  CANTILEVER  BEAMS,  NO  INERTIA  (V£L=CONST) 

Figure  A6-1: 
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:EMENT(INCHES)  WSPIACEMENT(INCHES)  DfSPLACEMENT(INCHES) 


TIME(SEC) 


TIME(SEC) 

VEL=  5.0  IN/SEC  ;  MASS=  1.0  ;  BEAM  LENGTH =  0.5  IN 


DYNAMIC  RESPONSE  OF  MESHING  CANTILEVER  BEAMS,  NO  INERTIA  (VEL=CONST) 


Figure  A6-2: 
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WSPLACEMENT(INCHES)  DISPLACEMENT(INCHES)  DISPLACEMENT(INCHES) 


.25 


.000  .005  .010  .015  .020 

TIME(SEC) 

VEL=  20.0  IN/SEC  ;  MASS=  1.0  ;  BEAM  LENGTH=  0.5  IN 

DYNAMIC  RESPONSE  OF  MESHING  CANTILEVER  BEAMS,  NO  INERTIA  (VEL-CONST) 

Figure  A6-4: 

218 


DISPIACEMENT(INCHES)  DISPLACEMENT(INCHES)  DfSPLACEMENT(INCHES) 


APPENDIX  7 

Dynamic  Response  Algorithm 
for  Meshing  Cantilever  Beam 

1.  Massless  Configuration 

2.  Inertia  of  Beam  Included 
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c*********************************************** 

C  THIS  PROGRAM  SOLVES  THE  EQUATION  OF  MOTION  * 

C  DESCRIBING  A  PAIR  OF  MESHING  CANTILEVER  * 

C  BEAMS  ATTACHED  TO  MOVEABLE  FOUNDATAION  * 

C  MASSES.  THE  SECOND  ORDER  DIFFERENTIAL  * 

C  EQUATION:  M(D2X/DT2)  +KX  -  P  IS  WRITTEN  * 

C  AS  A  FIRST  ORDER  DIFFEQ  AND  INTEGRATED  * 

C  USING  A  4TH  ORDER  "RUNGE-KUTTA"  INTEGRA-  * 

C  TION  ROUTINE.  * 

C*********************************************** 

EXTERNAL  FCT,OUTP 

DIMENSION  Y(4) , DERY(4) ,PRMT(5) ,AUX(8,4) 

REAL  P, L, BA, VEL, Ml , M2 , MASS , INERT, K1 , K2 , STIFF , MOD, STIFF1 
COMMON  P , L, VEL , MASS , INERT , K1 , K2 , STIFF , MOD , ZETA1 , ZETA2 , TF , DENOM 
WRITE (6,*)  'ENTER:  APPLIED  LOAD,  BEAM  LENGTH,  BASE  WIDTH' 
READ(5, *)  P,L, BA 

WRITE ( 6 , * )  'ENTER:  VELOCITY,  MASS#1,  MASS#2' 

READ ( 5 , * )  VEL, Ml, M2 
MOD-30. E6 

INERT* (1 ./12 . ) *BA*BA**3 
ZETA1-L/10. 

ZETA2=9.*L/10. 

K1- ( 3 . *MOD* INERT) / ( ZETA1**3 ) 

K2- ( 3 . *MOD* INERT) / ( ZETA2  *  *3 ) 


STIFF- (K1*K2 ) / (K1+K2 ) 

QTTPPI  aQ'PTPP 

WRITE ( 6 , * )  'APPLIED  LOAD  :',P 
WRITE ( 6 , * )  'BEAM  LENGTH  :',L 

WRITE ( 6 , * )  'BASE  WIDTH  :',BA 

WRITE ( 6 , * )  'MASS  #1  : • ,M1 

WRITE (6,*)  'MASS  #2  : • ,M2 

WRITE ( 6 , * )  'VELOCITY  : • ,VEL 


Y (1)— 0.0 

Y (2) -P/STIFF 

DERY (1)— 0. 5 

DERY ( 2 ) —0 . 5 

PRMT ( 1 ) =0 . 0 

PRMT ( 2 ) « ( . 8*L) /VEL 

PRMT ( 3 ) -PRMT ( 2 ) /5000 . 

WRITE (6,*)  'STARTING  TIME:  ',PRMT(1) 

WRITE ( 6 , * )  ! ENDING  TIME:  *,PRMT(2) 

WRITE ( 6 , * )  'TIME  INCREMENT:  ',PRMT(3) 

PRMT(4) -0.0001 
NDIM-2 

MASS- (M1*M2 ) / (M1+M2 ) 

TIME-0.0 

CALL  RKGS ( PRMT , Y , DERY , NDIM , IHLF , FCT , OUTP, AUX) 
******************************************************** 
******************************************************** 

SUBROUTINE  FCT (T, Y , DERY) 

DIMENSION  Y (4) , DERY ( 4 ) 

REAL  P, L, BA , VEL, Ml , M2 , MASS , INERT, K1 , K2 , STIFF, MOD 

COMMON  P , L , VEL, MASS , INERT , K1 , K2 , STI FF , MOD , ZETA1 , ZETA2 , TF , DENOM 

ZETA1- ( L/ 10 . ) +T*VEL 

ZETA2- ( 9 . *L/ 10 . ) -T*VEL 

Kl-(3. *MOD*INERT)/ (ZETA1**3) 

K2- ( 3 . *MOD* INERT) / ( ZETA2  *  *3) 

STIFF- (K1*K2 ) / (K1+K2 ) 

DERY ( 1 ) - ( P-STI FF* Y ( 2 ) ) /MASS 
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oonooonoooooooooooooooooo 


DERY(2)mY(l) 

DENOM- 1 . + ( ZETA1 * * 3/ ZETA2 * * 3 ) 

RETURN 

END 

****************************************************** 

****************************************************** 

SUBROUTINE  OUTP (T,  Y, DERY , IHLF , NDIM , PRMT) 

DIMENSION  Y(4) , DERY (4) 

COMMON  P, L, VEL, MASS , INERT , K1 , K2 , STIFF , MOD, ZETA1 , ZETA2 , TF, DENOM 

DELTA2-Y ( 2 ) /DENOM 

DELTAl-Y ( 2 ) -DELTA2 

WRITE (15,*)  T, Y ( 2 ) , DELTA1 , DELTA2 

RETURN 

END 

****************************************************** 

****************************************************** 

THE  SUBROUTINE  "RKGS"  (RUNGE-KUTTA)  IS  INCLUDED 
IN  THE  NEXT  PROGRAM 

RKG 

. RKG 

RKG 


SUBROUTINE  RKGS  RKG 

RKG 

PURPOSE  RKG 

TO  SOLVE  A  SYSTEM  OF  FIRST  ORDER  ORDINARY  DIFFERENTIAL  RKG 

EQUATIONS  WITH  GIVEN  INITIAL  VALUES.  RKG 

RKG 

USAGE  RKG 

CALL  RKGS  (PRMT, Y, DERY, NDIM, IHLF, FCT, OUTP, AUX)  RKG 

PARAMETERS  FCT  AND  OUTP  REQUIRE  AN  EXTERNAL  STATEMENT.  RKG 

RKG 

DESCRIPTION  OF  PARAMETERS  RKG 


PRMT  -  AN  INPUT  AND  OUTPUT  VECTOR  WITH  DIMENSION  GREATER  RKG 
OR  EQUAL  TO  5,  WHICH  SPECIFIES  THE  PARAMETERS  OF  RKG 
THE  INTERVAL  AND  OF  ACCURACY  AND  WHICH  SERVES  FOR  RKG 
COMMUNICATION  BETWEEN  OUTPUT  SUBROUTINE  (FURNISHED  RKG 


BY  THE  USER)  AND  SUBROUTINE  RKGS.  EXCEPT  PRMT (5)  RKG 
THE  COMPONENTS  ARE  NOT  DESTROYED  BY  SUBROUTINE  RKG 
RKGS  AND  THEY  ARE  RKG 

PRMT ( 1 ) -  LOWER  BOUND  OF  THE  INTERVAL  (INPUT),  RKG 

PRMT (2)-  UPPER  BOUND  OF  THE  INTERVAL  (INPUT),  RKG 

100:>TERVAL  (INPUT),  RKGS  220 

C  PRMT ( 2 ) -  UPPER  BOUND  OF  THE  INTERVAL 


t 
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C********************************************** 

C  THIS  PROGRAM  SOLVES  A  SYSTEM  OF  5  LINEAR  * 

C  DIFFERENTIAL  EQUATIONS  DESCRIBING  THE  * 

C  DYNAMIC  RESPONSE  OF  TWO  COUPLED  CANTILEVER  * 

C  BEAMS  WITH  A  MOVING  LOAD  BETWEEN  THEM.  * 

C ********************************************** 

c 

EXTERNAL  FCT,OUTP 

DIMENSION  FUNC(IO) ,  Y  (8 ) , DERY (8 ) , PRMT (5) ,AUX(8,8) ,A(5,5) 

REAL  MOD , INERT , Ml , M2 , MBl , MB2 , SI , S2 , OMEG , IPHI1 , IPHI2 
REAL  I2PHI1,I2PHI2, LAMBDA, LENG 

COMMON  IPHI1 , I PHI 2 , BETA, CONST, SI , S2 , Ml , M2 , MBl , MB2 , VEL 
$ , LAMBDA, ZETA1 , ZETA2 , I2PHI1 , I2PHI2 , LENG , PHI1 , PHI2 , DZZ 
WRITE (6,*)  'ENTER  APPLIED  LOAD ( lbs) , AND  SPEED  OF  MOVING  LOAD' 
READ(5, *)  P, VEL 

WRITE (6,*)  'ENTER  BEAMLENGTH ( In) , AND  BASE  THICKNESS (in) • 
READ(5, *)  LENG, THICK 

WRITE(6, *)  'ENTER  MASSES:  Ml, M2,  AND  DENSITY  OF  BEAM ( lbs/ in3 ) • 
READ(5, *)  Ml, M2, DENS 
MOD-30. E6 
LOO 

INERT* ( 1 ./12 . ) *THICK**4 
DENS- ( DENS/ 3 8 6 . 4 )  *THICK* * 2 
WRITE (6,*)  'APPLIED  LOAD:  ',P 
WRITE (6,*)  'VELOCITY  :  ',VEL 
WRITE (6, *)  'BEAM  LENGTH  :  ' , LENG 

WRITE (6,*)  'THICKNESS  :  ', THICK 

WRITE ( 6 , * )  'MASS/UNIT  :  ' , DENS 

C 

C****  ONE  TERM  OF  THE  NATURAL  MODE  SERIES  IS  USED  ****** 

C 

OMEG- ( 1 . 875**2) * (MOD*INERT/ (DENS*LENG**4) ) **0. 5 
BETA- ( ( (OMEG**2 ) *DENS) / (MOD* INERT) ) **0 . 25 
WRITE ( 6 , * )  'OMEGA-' ,OMEG, '  BETA-' , BETA 
C 

C****  EVALUATE  THE  CONSTANT  TERMS  IN  MODE  EQUATION  ****** 

C 

X- BETA* LENG 

51- SIN (X)-SINH(X) 

52- COS (X) +COSH(X) 

WRITE ( 6 , * )  'SI-', SI,'  S2— ' ,S2 
C 

C****  EVALUATE  CONSTANT  OF  INTEGRATION  OF  MODE  EQUATION  **** 

C****  INTEGRATE  USING  SIMPSON'S  RULE  -  100  INCREMENTS***** 

C 

XINC— LENG/100 . 

Z ETA-0.0 
XFACT-0.0 
DO  10  1-2,100,2 
DO  20  J-1,3 

ZETA-XINC* (1-2) +XINC* (J-l) 

X*  BETA  *  Z  ETA 

FUNC(J) -Sl**2* (SIN(X) **2+SINH(X) **2-2 . *SIN(X) *SINH(X) ) 
$+S2**2  * (COS (X) **2+COSH (X) **2-2 . *COS (X) *COSH (X) ) 

$+2. *S1*S2*(SIN(X) *COS(X) -SIN(X) *COSH(X) -COS(X) *SINH(X) 
$+SINH(X) *COSH(X) ) 

20  CONTINUE 

XFACT-XFACT+ (XINC/3 . ) * (FUNC ( 1) +4 . *FUNC (2 ) +FUNC ( 3 ) ) 

10  CONTINUE 
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CONST—  (l./(DENS*XFACT))**0.5 
WRITE (6, *)  'CONST* 1 , CONST 
C 

C****  EVALUATE  INTEGRAL  OF  MODE  SHAPE:  IPHI1-IPHI2**** 

C 

XINOLENG/IOO. 

ZETA-0.0 
XFACT-0 . 0 
DO  30  1-2,100,2 
DO  40  J-1,3 

ZETA-XINC* (1-2) +XINC* (J-l) 

X*Z  ETA*  BETA 

FUNC ( J) — (SI* (SIN (X) -SINH(X) )+S2*(COS(X) -COSH(X) ) ) 

40  CONTINUE 

XFACT-XFACT+ (XINC/3 . ) * (FUNC ( 1) +4 . *FUNC (2 ) +FUNC ( 3 ) ) 

30  CONTINUE 

IPHI1— DENS*XFACT*CONST 
IPHI2— IPHI1 

WRITE ( 6 , * )  'IPHI1— IPHI2— 1 , IPHI1 
C 

C****  EVALUATE  INTEGRAL  OF  (D2PHI/DZETA2) **2 :  I2PHI  **** 

C 

XINOLENG/IOO. 

ZETA-0.0 
XFACT-0. 0 
DO  50  1-2,100,2 
DO  60  J-1,3 

ZETA-XINC* ( 1-2 ) +XINC* (J-l) 

X-Z ETA* BETA 

FUNC (J)-(S1**2* (SIN (X) **2+SINH(X) **2 
$+2 . *SIN (X) *SINH (X) ) +S2**2* (COS (X) **2+COSH(X) **2 
$+2 . *COS (X) *COSH (X) ) +2 . *S1*S2* (SIN (X) *COS (X) +SIN (X) *COSH (X) 
$+COS (X) *SINH(X) +COSH(X) *SINH(X) ) ) 

60  CONTINUE 

XFACT— XFACT+ (XINC/3 . ) * (FUNC ( 1) +4 . *FUNC (2 ) +FUNC (3 ) ) 

50  CONTINUE 

1 2  PHI 1— XFACT*MOD* INERT *CONST** 2  *  BETA*  *4 
I2PHI2~I2PHI1 

WRITE (6, *)  'I2PHI1— I2PHI2— ' , I2PHI1 
C 

MB1-DENS*LENG 

MB2-MB1 

ZETA1-LENG/10. 

ZETA2— . 9  *LENG 
X1-ZETA1 *BETA 
X2-ZETA2 *  BETA 

PHI 1— CONST* (SI* (SIN(Xl) -SINH(Xl) )+S2*(C0S (XI) -COSH(Xl) ) ) 

PHI 2-CONST* (SI* (SIN (X2 ) -SINH (X2 ) ) +S2* (COS (X2 ) -COSH (X2) ) ) 

WRITE ( 6 , * )  'PHIl— ' , PHI1, '  PHI2-* , PHI2 

DZ-VEL 

DZZ-0.0 

D1PHI1-C0NST*BETA*DZ* (SI* (COS (XI) -COSH (XI) ) 

$-S2* (SIN (XI) +SINH(X1) ) ) 

D1PHI2— CONST*BETA* (-DZ) * (SI* (COS (X2) -COSH(X2) ) 

$-S2* (SIN (X2) +SINH(X2) ) ) 

C 

C****  DEFINE  THE  INITIAL  VALUES  OF  THE  DEPENDENT  VARIABLES  ** 

C 

Y(3)-P*PHI1/I2PHI1 
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Y (4 ) -P*PHX2/I2PHI2 
Y(l)-0.0 

Y(2)-¥(1)+PHI1*Y(3)+PHI2*Y(4) 

Y(7)-0.0 

Y(8)-0.0 

Y(5)-0.0 

Y (6) — Y (5)+Y (3) *D1PHI1+Y (4) *D1PHX2 
WRITE (6, *)  'Y(1)-',Y(1) 

WRITE (6,*)  *Y(2)-',Y(2) 

WRITE (6, *)  'Y(3)-',Y(3) 

WRITE (6,*)  *Y(4)-',Y(4) 

00  55  1-1,8 
DERY(X)— 1./8. 

55  CONTINUE 

PSMT(l)— 0.0 
PRMT(2)-.8*LENG/(VEL) 

PRMT ( 3 ) - ( PRMT ( 2 ) -PRMT ( 1 ) ) / 100 . 

PRMT (4) -0.0001 
NDIM-8 

CALL  RKGS ( PRMT , Y , DERY , NDIM , IHLF , FCT , OUTP , AUX ) 


SUBROUTINE  FCT(T, Y,DBRY) 

DIMENSION  Y(8) , DERY (8),A(5,5),B(5) 

REAL  MOD, INERT,M1,M2,MB1,MB2 ,S1,S2 ,OMEG, IPHI1, IPHI2 
REAL  I2PHI1 , I2PHI2 , LAMBDA, LENG 

COMMON  I PHI 1 , I PHI 2 , BETA , CONST , S 1 , S2 , Ml , M2 , MB1 , MB2 , VEL 
$ , LAMBDA , ZETA1 , ZETA2 , I2PHI1 , I2PHI2 , LENG, PHI1 , PHI2 , DZZ 
ZETA1— ( LENG/ 10 . ) +VEL*T 
ZETA2- ( . 9* LENG) -VEL*T 
XI— ZETA1*BETA 
X2-ZETA2  *BETA 
C 

C***«  EVALUATE  MODE  SHAPE  EQUATIONS  AT  EACH  LOAD  POSITION  *** 
C 

PHI1— CONST* (SI* (SIN (XI) -SINH(Xl) ) +S2* (COS(Xl) -COSH (XI) ) ) 
PHI 2 -CONST* (SI* (SIN(X2) -SINH(X2) ) +S2* (COS (X2) -COSH(X2) ) ) 

C 

IF(LC.GT.IO)  GO  TO  71 
WRITE(6, *)  '**********TIME-',T 
C****  EVALUATE  DERIVATIVES  OF  MODE  SHAPE  EQUATIONS  **** 

C 

71  DZ-VEL 
DZZ— 0 . 0 

D1PHI 1 -CONST*  BETA*  DZ * (SI* (COS (XI) -COSH (XI) ) 

$-S2* (SIN (XI) +SINH (XI) ) ) 

D1PHI 2 -CONST*  BETA* ( -DZ ) * ( S 1 * ( COS ( X2 ) -COSH ( X2 ) ) 

$-S2* (SIN (X2 ) +SINH (X2 ) ) ) 

D2PHI1-C0NST*BETA*(S1*(-BETA*DZ**2*SIN(X1)+DZZ*C0S(X1) 
$-BETA*DZ**2*SINH(Xl) -DZZ*C0SH(X1) ) -S2*(BETA*DZ**2*C0S(X1) 
$+DZZ*SIN(Xl)^BETA*DZ**2*COSH(Xl)+DZZ*SINH(Xl) ) ) 

D2 PHI 2— CONST* BETA* (Sl*(-BETA*DZ**2*SIN(X2)+DZZ*COS (X2) 
$-BETA*DZ**2*SINH(X2) -DZZ*COSH(X2) ) -S2*(BETA*DZ**2*COS(X2) 
$>DZZ*SIN(X2)+BETA*DZ**2*COSH(X2)+DZZ*SINH(X2) ) ) 

C 

C****  DEFINE  MATRIX  COEFFICIENTS  A[5,5]  **** 

C 


A(l, 1)-M1+MB1 
A(2, 1)— 0.0 


225 


A(3,1)-IPHI1 
A(4,l)-0.0 
A<5,1)  —  1. 

A(l,2)-0.0 
A(2,2)-M2+MB2 
A(3,2)*0.0 
A(4,2)— IPHI2 
A(5,2)-1.0 
A(1,3)*IPHI1 
A(2,3)-0.0 
A(3,3)*1.0 
A(4,3)«0.0 
A(5f3)— PHI1 
A(l,4)*0.0 
A(2,4)«- IPHI2 
A(3,4)-0.0 
A<4,4)-1.0 
A(5,4)»- PHI2 
A(l,5)»1.0 
A(2,5)«- 1.0 
A(3,5)-PHI1 
A(4,5)-PHI2 
A(5,5)-0.0 
C 

C****  DEFINE  RIGHT  HAND  SIDE  B[5]  ******* 

C 

B(l)— P  . 

B(2)-P 

B(3)«- I2PHI1*Y(3) 

B(4)«-  I2PHI2*Y(4) 

B(5) -D2PHI1*Y (3) +D2PHI2*Y (4) +2 . *D1PHI1*Y (7) +2 . *D1PHI2*Y (8) 
IF(LC.GT.IO)  GO  TO  1 
DO  63  1-1,5 

63  WRITE(6,*)  A(I,1),A(I,2),A(I,3),A(I,4),A(I,5),B(I) 

1  LC-LC+1 

N— 5 

CALL  SIMQ(A,B,N,KS) 

IF(KS.EQ.l)  WRITE  (6,  *)  'NO  SOLUTION!  Ml  I!  • 

IF(LC.GT.IO)  GO  TO  64 
WRITE(6, *)  'BACK  FROM  SIMQ  [B] ' 

DO  64  1-1,5 
WRITE (6, *)  B(I) 

64  CONTINUE 
DERY (1) — Y(5) 

DERY(2)-Y(6) 

DERY (3) -Y (7) 

DERY ( 4 ) — Y ( 8 ) 

DERY(5)-B(1) 

DERY  (6)  — B(2) 

DERY (7)— B(3) 

DERY (8)— B(4) 

LAMBDA-B(5) 

IF(LC.GT.IO)  GO  TO  62 
DO  61  1-1,8 

61  WRITE(6, *)  'Y(',I, »)-',Y(I) 

62  RETURN 
END 
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SUBROUTINE  OUTP(T, Y,DERY, IHLF,NDIM, PRMT) 

DIMENSION  Y ( 8 ) , DERY (  8 )  , PRMT ( 5 ) 

REAL  MOD, INERT , Ml , M2 , MB1 , MB2 , SI , S2 , OMEG, IPHI1 , IPKI2 
REAL  I 2 PHI 1,1 2 PHI 2, LAMBDA, LENG 

COMMON  I PHI 1 , I PHI 2 , BETA , CONST , S 1 , S2 , Ml , M2 , MB1 , MB2 , VEL 
$ , LAMBDA, ZETA1 , ZETA2 , I2PHI1 , I2PHI2 , LENG, PHI1 , PHI2 , DZZ 
ZETA1- (LENG/10 . ) +VEL*T 
ZETA2- ( . 9*LENG) -VEL*T 
X1=ZETA1*BETA 
X2=ZETA2 *BETA 

PHI l»CONST * (SI* (SIN (XI) -SINH(Xl) ) +S2* (COS (XI) -COSH (XI) ) ) 

PHI2“CONST* (SI* (SIN (X2) -SINH (X2 ) ) +S2* (COS (X2 ) -COSH(X2) ) ) 

WRITE ( 15, *)  T,Y(2)-Y(1) ,Y(3) *PHI1, Y (4) *PHI2 

WRITE ( 6 , * )  ‘ACCEL*' ,DZZ 

RETURN 

END 

************************************************ 

************************************************ 


SUBROUTINE  SIMQ 
PURPOSE 

OBTAIN  SOLUTION  OF  A  SET  OF  SIMULTANEOUS  LINEAR  EQUATIONS, 
AX«B 

USAGE 

CALL  SIMQ(A,B,N,KS) 

DESCRIPTION  OF  PARAMETERS 

A  -  MATRIX  OF  COEFFICIENTS  STORED  COLUMNWISE.  THESE  ARE 
DESTROYED  IN  THE  COMPUTATION.  THE  SIZE  OF  MATRIX  A  IS 
N  BY  N. 

B  -  VECTOR  OF  ORIGINAL  CONSTANTS  (LENGTH  N) .  THESE  ARE 
REPLACED  BY  FINAL  SOLUTION  VALUES,  VECTOR  X. 

N  -  NUMBER  OF  EQUATIONS  AND  VARIABLES.  N  MUST  BE  .GT.  ONE. 
KS  -  OUTPUT  DIGIT 

0  FOR  A  NORMAL  SOLUTION 
1  FOR  A  SINGULAR  SET  OF  EQUATIONS 


REMARKS 

MATRIX  A  MUST  BE  GENERAL. 

IF  MATRIX  IS  SINGULAR  ,  SOLUTION  VALUES  ARE  MEANINGLESS. 

AN  ALTERNATIVE  SOLUTION  MAY  BE  OBTAINED  BY  USING  MATRIX 
INVERSION  (MINV)  AND  MATRIX  PRODUCT  (GMPRD) . 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

METHOD 

METHOD  OF  SOLUTION  IS  BY  ELIMINATION  USING  LARGEST  PIVOTAL 
DIVISOR.  EACH  STAGE  OF  ELIMINATION  CONSISTS  OF  INTERCHANGING 
ROWS  WHEN  NECESSARY  TO  AVOID  DIVISION  BY  ZERO  OR  SMALL 
ELEMENTS. 

THE  FORWARD  SOLUTION  TO  OBTAIN  VARIABLE  N  IS  DONE  IN 
N  STAGES.  THE  BACK  SOLUTION  FOR  THE  OTHER  VARIABLES  IS 
CALCULATED  BY  SUCCESSIVE  SUBSTITUTIONS.  FINAL  SOLUTION 
VALUES  ARE  DEVELOPED  IN  VECTOR  B,  WITH  VARIABLE  1  IN  B(l), 
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C  VARIABLE  2  IN  B(2) , . .  VARIABLE  N  IN  B(N)  . 

C  IF  NO  PIVOT  CAN  BE  FOUND  EXCEEDING  A  TOLERANCE  OF  0.0, 

C  THE  MATRIX  IS  CONSIDERED  SINGULAR  AND  KS  IS  SET  TO  1.  THIS 
C  TOLERANCE  CAN  BE  MODIFIED  BY  REPLACING  THE  FIRST  STATEMENT. 
C 

C  . 

c 

SUBROUTINE  SIMQ(A,B,N,KS) 

DIMENSION  A(l) ,  B(l) 

C 

C  FORWARD  SOLUTION 

C 

TOL-O.O 
KS-0 
JJ— N 

DO  65  J-1,N 
JY-J+1 
JJ-JJ+N+1 
BIGA-0 
IT-JJ-J 
DO  30  I-J,N 
C 

C  SEARCH  FOR  MAXIMUM  COEFFICIENT  IN  COLUMN 

C 

IJ-IT+I 

IF (ABS (BIGA) -ABS (A(IJ) ) )  20,30,30 
20  BIGA-A(IJ) 

I MAX* I 
30  CONTINUE 
C 

C  TEST  FOR  PIVOT  LESS  THAN  TOLERANCE  (SINGULAR  MATRIX) 

C 

IF (ABS ( BIGA) -TOL)  35,35,40 
35  KS-1 
RETURN 
C 

C  INTERCHANGE  ROWS  IF  NECESSARY 

C 

40  Il*J+N*(J-2) 

IT-IMAX-J 
DO  50  K«J,N 

11- Il+N 

12- I1+IT 
SAVE-A(Il) 

A(I1)-A(I2) 

A (12) -SAVE 

C 

C  DIVIDE  EQUATION  BY  LEADING  COEFFICIENT 

C 

50  A(I1) — A(I1) /BIGA 
SAVE— B ( IMAX) 

B(IMAX)— B(J) 

B(J) -SAVE/ BIGA 
C 

C  ELIMINATE  NEXT  VARIABLE 

C 

IF(J-N)  55,70,55 
55  IQS— N* (J-l) 

DO  65  IX— JY , N 
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IXJ-IQS+IX 

IT-J-IX 

DO  60  JX-JY,N 

IXJX-N*(JX-1)+IX 

JJX-IXJX+IT 

60  A(IXJX)-A(IXJX)-(A(IXJ)*A(JJX)) 
65  B(IX)«B(IX)-(B(J)*A(IXJ)) 

BACK  SOLUTION 

70  NY-N-1 

it-n*n 

DO  80  J-1,NY 
IA-IT-J 
IB»N-J 
IC-N 

DO  80  K«1,J 

B(IB) -B(IB) -A(IA) *B(IC) 

IA-IA-N 
80  IC-IC-1 
RETURN 
END 


SUBROUTINE  RKGS 

PURPOSE 

TO  SOLVE  A  SYSTEM  OF  FIRST  ORDER  ORDINARY  DIFFERENTIAL 

EQUATIONS  WITH  GIVEN  INITIAL  VALUES. 

USAGE 

CALL  RKGS  ( PRMT , Y , DERY , NDIM , IHLF , FCT , OUTP , AUX ) 

PARAMETERS  FCT  AND  OUTP  REQUIRE  AN  EXTERNAL  STATEMENT. 

DESCRIPTION  OF  PARAMETERS 

PRMT  -  AN  INPUT  AND  OUTPUT  VECTOR  WITH  DIMENSION  GREATER 
OR  EQUAL  TO  5,  WHICH  SPECIFIES  THE  PARAMETERS  OF 
THE  INTERVAL  AND  OF  ACCURACY  AND  WHICH  SERVES  FOR 
COMMUNICATION  BETWEEN  OUTPUT  SUBROUTINE  (FURNISHED 
BY  THE  USER)  AND  SUBROUTINE  RKGS.  EXCEPT  PRMT (5) 

THE  COMPONENTS  ARE  NOT  DESTROYED  BY  SUBROUTINE 
RKGS  AND  THEY  ARE 

PRMT(l) -  LOWER  BOUND  OF  THE  INTERVAL  (INPUT), 

PRMT (2)-  UPPER  BOUND  OF  THE  INTERVAL  (INPUT), 

PRMT (3)-  INITIAL  INCREMENT  OF  THE  INDEPENDENT  VARIABLE 
(INPUT) , 

PRMT (4)-  UPPER  ERROR  BOUND  (INPUT).  IF  ABSOLUTE  ERROR  IS 
GREATER  THAN  PRMT (4) ,  INCREMENT  GETS  HALVED. 

IF  INCREMENT  IS  LESS  THAN  PRMT (3)  AND  ABSOLUTE 
ERROR  LESS  THAN  PRMT(4)/50,  INCREMENT  GETS  DOUBLED. 
THE  USER  MAY  CHANGE  PRMT (4}  BY  MEANS  OF  HIS 
OUTPUT  SUBROUTINE. 

PRMT (5)-  NO  INPUT  PARAMETER.  SUBROUTINE  RKGS  INITIALIZES 
PRMT(5)«0.  IF  THE  USER  WANTS  TO  TERMINATE 
SUBROUTINE  RKGS  AT  ANY  OUTPUT  POINT,  HE  HAS  TO 
CHANGE  PRMT (5)  TO  NON-ZERO  BY  MEANS  OF  SUBROUTINE 
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OUTP.  FURTHER  COMPONENTS  OF  VECTOR  PRMT  ARE 
FEASIBLE  IF  ITS  DIMENSION  IS  DEFINED  GREATER 
THAN  5.  HOWEVER  SUBROUTINE  RKGS  DOES  NOT  REQUIRE 
AND  CHANGE  THEM.  NEVERTHELESS  THEY  MAY  BE  USEFUL 
FOR  HANDING  RESULT  VALUES  TO  THE  MAIN  PROGRAM 
(CALLING  RKGS)  WHICH  ARE  OBTAINED  BY  SPECIAL 
MANIPULATIONS  WITH  OUTPUT  DATA  IN  SUBROUTINE  OUTP. 

Y  -  INPUT  VECTOR  OF  INITIAL  VALUES.  (DESTROYED) 

LATERON  Y  IS  THE  RESULTING  VECTOR  OF  DEPENDENT 
VARIABLES  COMPUTED  AT  INTERMEDIATE  POINTS  X. 

DERY  -  INPUT  VECTOR  OF  ERROR  WEIGHTS.  (DESTROYED) 

THE  SUM  OF  ITS  COMPONENTS  MUST  BE  EQUAL  TO  1. 
LATERON  DERY  IS  THE  VECTOR  OF  DERIVATIVES,  WHICH 
BELONG  TO  FUNCTION  VALUES  Y  AT  A  POINT  X. 

NDIM  -  AN  INPUT  VALUE,  WHICH  SPECIFIES  THE  NUMBER  OF 
EQUATIONS  IN  THE  SYSTEM. 

IHLF  -  AN  OUTPUT  VALUE,  WHICH  SPECIFIES  THE  NUMBER  OF 

BISECTIONS  OF  THE  INITIAL  INCREMENT.  IF  IHLF  GETS 
GREATER  THAN  10,  SUBROUTINE  RKGS  RETURNS  WITH 
ERROR  MESSAGE  IHLF=11  INTO  MAIN  PROGRAM.  ERROR 
MESSAGE  IHLF=12  OR  IHLF=13  APPEARS  IN  CASE 
PRMT ( 3 ) =0  OR  IN  CASE  SIGN(PRMT(3) ) .NE.SIGN(PRMT(2) - 
PRMT ( 1 ) )  RESPECTIVELY. 

FCT  -  THE  NAME  OF  AN  EXTERNAL  SUBROUTINE  USED.  THIS 

SUBROUTINE  COMPUTES  THE  RIGHT  HAND  SIDES  DERY  OF 
THE  SYSTEM  TO  GIVEN  VALUES  X  AND  Y.  ITS  PARAMETER 
LIST  MUST  BE  X , Y , DERY .  SUBROUTINE  FCT  SHOULD 
NOT  DESTROY  X  AND  Y. 

OUTP  -  THE  NAME  OF  AN  EXTERNAL  OUTPUT  SUBROUTINE  USED. 

ITS  PARAMETER  LIST  MUST  BE  X, Y, DERY, IHLF, NDIM, PRMT. 
NONE  OF  THESE  PARAMETERS  (EXCEPT,  IF  NECESSARY, 

PRMT ( 4 ) , PRMT ( 5 ) , . . . )  SHOULD  BE  CHANGED  BY 
SUBROUTINE  OUTP.  IF  PRMT (5)  IS  CHANGED  TO  NON-ZERO, 
SUBROUTINE  RKGS  IS  TERMINATED. 

AUX  -  AN  AUXILIARY  STORAGE  ARRAY  WITH  8  ROWS  AND  NDIM 
COLUMNS . 

REMARKS 

THE  PROCEDURE  TERMINATES  AND  RETURNS  TO  CALLING  PROGRAM,  IF 

(1)  MORE  THAN  10  BISECTIONS  OF  THE  INITIAL  INCREMENT  ARE 
NECESSARY  TO  GET  SATISFACTORY  ACCURACY  (ERROR  MESSAGE 
IHLF=11) , 

(2)  INITIAL  INCREMENT  IS  EQUAL  TO  0  OR  HAS  WRONG  SIGN 
(ERROR  MESSAGES  IHLF=12  OR  IHLF=13) , 

(3)  THE  WHOLE  INTEGRATION  INTERVAL  IS  WORKED  THROUGH, 

(4)  SUBROUTINE  OUTP  HAS  CHANGED  PRMT (5)  TO  NON-ZERO. 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

THE  EXTERNAL  SUBROUTINES  FCT (X, Y , DERY)  AND 

OUTP (X,Y, DERY, IHLF, NDIM, PRMT)  MUST  BE  FURNISHED  BY  THE  USER. 

METHOD 

EVALUATION  IS  DONE  BY  MEANS  OF  FOURTH  ORDER  RUNGE-KUTTA 

FORMULAE  IN  THE  MODIFICATION  DUE  TO  GILL.  ACCURACY  IS 

TESTED  COMPARING  THE  RESULTS  OF  THE  PROCEDURE  WITH  SINGLE 

AND  DOUBLE  INCREMENT. 

SUBROUTINE  RKGS  AUTOMATICALLY  ADJUSTS  THE  INCREMENT  DURING 

THE  WHOLE  COMPUTATION  BY  HALVING  OR  DOUBLING.  IF  MORE  THAN 

10  BISECTIONS  OF  THE  INCREMENT  ARE  NECESSARY  TO  GET 
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SATISFACTORY  ACCURACY,  THE  SUBROUTINE  RETURNS  WITH 
ERROR  MESSAGE  IHLF-11  INTO  MAIN  PROGRAM. 

TO  GET  FULL  FLEXIBILITY  IN  OUTPUT,  AN  OUTPUT  SUBROUTINE 
MUST  BE  FURNISHED  BY  THE  USER. 

FOR  REFERENCE,  SEE 

RALSTON/WXLF,  MATHEMATICAL  METHODS  FOR  DIGITAL  COMPUTERS, 
WILEY,  NEW  YORK/LONDON,  1960,  PP. 110-120. 


SUBROUTINE  RJCGS  (PRMT,  Y,  DERY,NDIM,  IHLF,  FCT.OUTP,  AUX) 


DIMENSION  Y'1) .  DERY(l) , AUX (8, 1) ,A(4),B(4),C(4) ,PRMT(1) 
DO  1  I-1.ND1M 

1  AUX (8, I J-. 06666667* DERY (I) 

X-PRMT ( 1 ) 

XEND-PRMT ( 2 ) 

H-PRMT ( 3 ) 

PRMT(5)-0. 

CALL  FCT ( X , Y , DERY ) 

ERROR  TEST 

IF(H* (XEND-X) ) 38 , 37 , 2 

PREPARATIONS  FOR  RUNGE-KUTTA  METHOD 

2  A(l)-.5 
A(2)-. 2928932 
A(3)-l. 707107 
A(4)-. 1666667 
B(l)-2. 

B(2)-l. 

B(3)-l. 

B(4)-2. 

C(l)“ • 5 
C(2)«. 2928932 
C(3)«l. 707107 
C(4)=. 5 

PREPARATIONS  OF  FIRST  RUNGE-KUTTA  STEP 
DO  3  I-1,NDIM 
AUX (1,1) *Y ( I ) 

AUX (2,1) *DERY ( I ) 

AUX(3,I)«0. 

3  AUX (6,1) *0 . 

IREC-0 
H«H+H 
IHLF— 1 
ISTEP-0 
IEND-0 


START  OF  A  RUNGE-KUTTA  STEP 

4  IF ( (X+H-XEND) *H) 7,6,5 

5  H-XEND-X 

6  IEND-1 

RECORDING  OF  INITIAL  VALUES  OF  THIS  STEP 

7  CALL  OUTP(X,Y, DERY, IREC,NDIM, PRMT) 
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oo  oo  no  oo  oooo  ooo 


IF(PRMT(5) ) 40, 8 , 40 

8  ITEST-0 

9  ISTEP-ISTEP+1 


START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
J«1 

10  AJ=A(J) 

BJ=B ( J) 

CJ=C(J) 

DO  11  I-1,NDIM 
R1«H*DER¥(I) 

R2«AJ*(R1-BJ*AUX(6,I) ) 

¥ ( I ) — Y ( I ) +R2 

11  AUX(6,I)=AUX(6,I)+R2-CJ*R1 
IF(J-4) 12,15, 15 

12  J=J+1 
IF(J-3)13,14,13 

13  X«X+.5*H 

14  CALL  FCT(X,Y,DERY) 

GOTO  10 

END  OF  INNERMOST  RUNGE-KUTTA  LOOP 


TEST  OF  ACCURACY 
15  IF(ITEST) 16, 16,20 


IN  CASE  ITEST=0  THERE  IS  NO  POSSIBILITY  FOR  TESTING  OF  ACCURACY 

16  DO  17  I=1,NDIM 

17  AUX (4,1) = Y ( I ) 

ITEST-1 

ISTEP=ISTEP+ISTEP-2 

18  IHLF-IHLF+1 
X-X-H 
H«.5*H 

DO  19  I-1,NDIM 
Y ( I ) “AUX (1,1) 

DERY ( I ) "AUX (2,1) 

19  AUX (6, I) -AUX (3, I) 

GOTO  9 

IN  CASE  ITEST-1  TESTING  OF  ACCURACY  IS  POSSIBLE 

20  IMOD— ISTEP/2 

IF ( ISTEP-IMOD-IMOD) 21,23,21 

21  CALL  FCT(X, Y, DERY) 

DO  22  I— 1,NDIM 
AUX (5,1) -Y ( I ) 

22  AUX (7,1) —DERY ( I ) 

GOTO  9 


COMPUTATION  OF  TEST  VALUE  DELT 

23  DELT-0. 

DO  24  I— 1,NDIM 

24  DELT-DELT+AUX (8,1) *ABS (AUX ( 4 , I ) -Y ( I ) ) 
IF(DELT-PRMT(4) ) 28 , 28, 25 


ERROR  IS  TOO  GREAT 
25  IF(IHLF-IO) 26,36,36 
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26  DO  27  I-1,NDIM 

27  AUX(4,I)«AUX(5,I) 

ISTEP-ISTEP+ISTEP-4 

X-X-H 

I END-0 
GOTO  18 
C 

C  RESULT  VALUES  ARE  GOOD 

28  CALL  FCT(X,Y,DERY) 

DO  29  I— 1 ,NDIM 
AUX(1,I)»Y(I) 

AUX (2,1) — DERY ( I ) 

AUX (3,1) — AUX (6,1) 

Y ( I ) —AUX (5,1) 

29  DERY ( I ) —AUX (7,1) 

CALL  OUTP (X-H , Y , DERY , IHLF , NDIM , PRMT) 
IF(PRMT(5) ) 40, 30, 40 

30  DO  31  I— 1, NDIM 
Y ( I ) —AUX (1,1) 

31  DERY ( I ) —AUX (2,1) 

IREC-IHLF 

IF(IEND)32,32,39 

C 

C  INCREMENT  GETS  DOUBLED 

32  IHLF— IHLF- 1 
ISTEP—ISTEP/ 2 
H“H+H 

C  TO  STOP  AT  INPUT  DELT:  IF (IHLF)  4,33,33 

IF (IHLF)  4,33,33 

33  IM0D=*ISTEP/2 

IF ( ISTEP-IMOD-IMOD) 4,34,4 

34  IF(DELT-. 02*PRMT(4) ) 35 , 35 , 4 

35  IHLF— IHLF- 1 
ISTEP— ISTEP/2 
H-H+H 

GOTO  4 
C 
C 

C  RETURNS  TO  CALLING  PROGRAM 

36  IHLF— 11 

CALL  FCT(X, Y/JERY) 

GOTO  39 

37  IHLF— 12 
GOTO  39 

38  IHLF— 13 

39  CALL  OUTP(X,Y, DERY, IHLF, NDIM, PRMT) 

40  RETURN 
END 
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